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1.  Introduction 


The  capability  to  calculate  the  burning  rate  of  propellants  from  their  ingredients  has  long  been 
recognized  as  desirable  though,  until  recently,  not  realizable.  Even  today,  this  capability  is  limited 
to  classes  of  energetic  materials  that  have  been  fairly  extensively  studied  experimentally.  So  far, 
the  goal  of  computing  all  the  properties  of  propellants  with  ingredients  that  have  never  actually 
been  synthesized  is  still  beyond  our  reach.  However,  considerable  progress  towards  obtaining 
properties  of  notional  materials  such  as  heats  of  formation,  density,  and  detonation  sensitivity  has 
been  made  in  the  last  10  years.  Over  the  same  period  of  time,  fledgling  ability  to  compute  the 
burning  rate  for  unstudied  formulations  within  known  classes  of  propellants  has  emerged. 

What  is  it  that  makes  the  buming-rate  calculations  so  difficult?  The  combustion  of  energetic 
materials  involves  coupling  physical  and  chemical  phenomena  in  the  condensed  phase,  the  gas 
phase,  and  at  the  interface  between  the  two.  Our  greatest  knowledge  and  experience  is  centered  in 
the  gas  phase.  We  now  know  that  many  dozens  of  chemical  species  reacting  by  many  hundreds  of 
elementary  reactions  are  at  play  there  in  concert  with  the  physical  processes  of  molecular 
diffusion,  convection,  and  thermal  conduction.  Even  by  itself,  the  gas  phase  presents  us  with  a 
daunting  array  of  nonlinear  phenomena  to  understand.  Fortunately,  general  scientific  progress 
over  the  last  several  decades  has  armed  us  with  the  experimental  and  theoretical  tools  to  approach 
this  problem  in  a  systematic  way.  Though  many  uncertainties  remain,  the  conceptual  means,  if  not 
always  the  resources,  are  available  to  resolve  them.  Such  is  not  the  case  in  dealing  with  the 
condensed-phase  and  interfacial  processes.  In  these  cases,  we  truly  do  not  know  what  we  do  not 
know. 

It  would  seem  that  the  best  hope  for  breaking  this  impasse  lies  in  molecular  dynamics  (MD) 
simulation  using  reactive  potentials,  but  this  approach  is  fairly  said  to  be  in  its  infancy.  Thus  far, 
only  continuum-mechanics  models  have  been  developed  to  describe  the  2-  or  3-phase  combustion 
process,  and  arguments  are  easily  mounted  to  suggest  that  one  will  always  want  to  use  a  continuum 
model  to  describe  the  gas  phase.  Yet,  I  believe  that  mating  these  two  descriptions  will  not  be  a 
trivial  task.  Anticipation  of  the  coming  need  for  communication  between  practitioners  of  these 
two  approaches  is  what  motivates  this  discourse.  Accordingly,  emphasis  is  placed  on  the 
mathematical  framework  and  concepts  behind  the  continuum-mechanics  approach,  the  history  of 
model  development,  their  current  status,  and  unsolved  problems. 
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2.  Phenomena 


Before  approaching  the  issues  surrounding  the  mathematical  modeling  of  propellant  burning  rates, 
one  should  become  familiar  with  the  range  of  phenomenological  behaviors  associated  with  the 
combustion  of  propellants  and  propellant  ingredients.  Solid  propellants  are  of  two  basic  types: 
homogeneous  and  composite.  Single-  and  double-base  propellants  are  examples  of  homogeneous 
propellants.  Homogeneous  propellants  are  considered  to  be  well  mixed  on  a  molecular  scale, 
although  this  may  not  be  strictly  true  because  of  the  practical  constraints  on  mixing.  Single -base 
propellant  consists  almost  exclusively  of  nitrocellulose  (NC),  which  becomes  increasingly  less 
soluble  in  the  manufacturing  solvents  as  its  nitration  level  increases  above  ~13%N,  leading  to  a 
very  viscous  liquid  during  mixing  and  extrusion.  My  first  experience  in  measuring  burning  rates 
was  with  6-in-long  strands  of  M10.  Based  on  the  measurement  technique  and  degree  of  control 
over  known  variables  such  as  pressure,  I  was  expecting  a  standard  error  of  about  1%.  Instead,  the 
standard  deviation  in  measured  burning  rates  using  dozens  of  1 -in-long  specimens  was  as  large  as 
25%.  These  large  deviations  arose  from  inhomogeneities  in  the  propellant  material  itself  due  to 
imperfect  mixing  of  the  highly  viscous  feedstock  material.  This  was  a  most  inauspicious 
introduction  to  my  new  field!  Fortunately,  this  would  become  the  worst  case  I  would  ever 
encounter  in  my  career.  Most  other  homogeneous  propellants  benefit  from  the  plasticizing 
properties  of  nitroglycerin  (NG)  and  burn  very  reproducibly.  However,  one  always  must  be  aware, 
particularly  with  experimental  propellants,  that  manufacturing  procedures  may  influence  the 
combustion  properties  in  unexpected  and  nonreproducible  ways. 

Composite  propellants  comprise  a  heterogeneous  mixture  of  crystalline  oxidizers  and 
polymeric  binders.  Double -base  propellant  with  crystalline  nitroguanidine  (NQ)  added 
(M30)  is  an  example.  Only  in  the  last  15  years  have  composite  propellants  involving 
cyclotrimethylenetrinitramine  (RDX)  come  into  use  for  guns.  The  use  of  such  a  secondary 
explosive  as  a  major  ingredient  in  gun  propellant  required  considerable  testing  to  pass  safety 
criteria  for  use  in  the  field.  Composites  utilizing  RDX,  cyclotetramethylenetetranitramine  (HMX), 
and  2,4,6,8,10,12-hexanitrohexaazaisowurtzitane  (CL20)  coupled  with  various  energetic  polymeric 
binders  are  currently  under  active  investigation. 

The  propellant  type  can  have  important  consequences  for  modeling  the  burning  rate.  In  many 
cases,  it  is  found  that  the  oxidizer  particle  size  strongly  influences  the  burning  rate.  If  this  is  so,  a 
one-dimensional  (1-D)  model  may  not  be  suitable.  On  the  other  hand,  observations  of  burning 
surfaces  of  nitramine-composite  propellants  ( 1 )  at  pressures  on  the  order  of  1  MPa  indicate  that  a 
melt  layer  exists  at  the  surface.  This  melt  layer  provides  an  opportunity  for  the  solid  ingredients  to 


Propellant  formulations  are  given  in  the  Appendix. 
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become  intimately  mixed  prior  to  gasification,  potentially  restoring  a  1-D  character  to  the 
combustion  phenomena. 

Despite  their  uniform  morphology,  one  cannot  necessarily  assume  that  homogeneous  propellants 
bum  one-dimensionally.  The  following  description  of  double-base  propellant  combustion  is 
enough  to  give  pause  to  the  most  intrepid  model  builder:  “It  can  be  observed  visually  that  the 
burning  surface  exhibits  a  wave-like  mode  of  consumption.  It  appears  as  if  glowing  filaments  of 
carbonaceous  material  periodically  move  over  the  surface  consuming  a  thin  layer  of  propellant. 
Cine  photography  of  the  propellant  surface  shows  that  a  smooth  area  on  this  surface  appears  to 
darken  and  to  roughen.  This  area  is  then  consumed  by  a  wave  of  combustion  which  moves  across 
the  surface  and  leaves  behind  a  network  of  carbon  filaments  which  are  blown  off  the  surface  by  the 
steady  evolution  of  gas.  ...  The  consumption  of  the  reacting  surface  layer  leaves  a  smooth  surface 
which  subsequently  repeats  this  sequence.  It  seems  probable  that  the  overall,  average  rate  of 
burning  is  the  result  of  two  components:  a  steady  rate  analogous  to  the  burning  of  a  liquid 
propellant  and  a  surface  or  condensed  phase  reactive  wave  moving  laterally  across  the  surface” 

(2).  It  is  clear  that  any  1-D  model  of  homogeneous  or  composite  propellant  combustion  is  to  be 
understood  as  an  idealization  of  the  phenomena  and  not  as  an  exact  description. 

Despite  all  of  the  aforementioned  complexities,  propellants  burn,  at  least  on  a  macroscopic  scale, 
in  “parallel  layers”  (i.e.,  in  a  direction  normal  to  the  local  surface  curvature).  This  is  the  property 
that  allows  for  interior-ballistic  control  of  the  net  gasification  rate  through  intricate  and  ingenious 
propellant-grain  geometries.  It  is  also  the  property  that  gives  the  model  builder  some  basis  for 
hope.  Another  source  of  encouragement  is  that  most  propellants  burn  with  a  very  simple  power- 
law  pressure  dependence,  despite  radical  changes  in  the  appearance  of  the  visible  flame  attached  to 
the  burning  surface.  At  pressures  below  ~1  MPa,  there  is  no  visible  flame  above  the  burning 
surface  of  homogeneous  and  many  composite  propellants.  As  the  pressure  increases,  a  weak  flame 
appears  >1  cm  distant  from  the  surface.  With  further  pressure  increases,  the  visible  flame  resides 
closer  and  closer  to  the  surface,  under  steady-state,  constant-pressure  conditions.  At  >10  MPa,  the 
flame  appears  attached  directly  to  the  surface.  In  most  cases,  this  evolving  flame  behavior  does 
not  perturb  the  power-law  pressure  dependence  of  the  burning  rate.  The  nonluminous  zone 
between  the  visible  flame  and  the  surface  is  known  as  the  dark  zone,  which  is  now  known  to  be  a 
consequence  of  the  relatively  slow  reduction  of  NO  to  N?  (and  slow  reactions  involving  HCN  in 
the  case  of  nitramines).  All  this  suggests  that  buming-rate  control  resides  in  a  thin  gas  layer  close 
to  the  surface,  an  inference  that  further  encourages  the  1-D  idealization. 


3.  Concepts 


To  date,  there  have  been  no  MD  models  of  multiphase  combustion.  This  may  change  in  the 
distant  future;  however,  a  nearer-term  prospect  is  that  MD  submodels  of  the  condensed-phase 
and/or  interface  processes  will  be  developed  and  that  these  submodels  will  be  used  to  supplement 
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the  continuum-mechanics  combustion  models.  Anticipating  the  need  to  merge  these  two 
viewpoints,  we  develop  in  this  section  the  mathematical  formulation  of  the  continuum  paradigm. 

Studying  Figure  1  will  provide  the  reader  with  some  appreciation  of  the  range  of  phenomena  that 
are  involved  in  the  combustion  of  propellants.  No  model  exists  which  has  treated  all  of  the 
processes  suggested  in  the  figure;  the  list  is  intended  as  a  conceptual  transition  from  nature  to 
mathematics.  Some  propellants  are  known  to  exhibit  a  liquid  layer  at  the  burning  surface  and 
some  may  not.  I  shall  assume  that,  in  the  general  case,  a  3-phase  (solid,  liquid,  gas)  problem  must 
be  solved.  The  photograph  in  Figure  1  is  of  an  RDX  composite  propellant  (M43)  burning  in  the 
steady  state  at  a  pressure  of  1.6  MPa.  In  this  case,  the  dark  zone  is  clearly  defined,  and  one  can  see 
that  the  idealization  of  a  1-D  semi-infinite  solid  does  not  appear  to  be  unreasonable. 


GAS-PHASE  FLAME 

♦  elementary  reactions 

♦  thermal  conduction 

♦  convection 

♦  molecular  diffusion 

♦  multi-component  transport 

♦  thermal  diffusion 


(  MULTI-COMPONENT 
EVAPORATION  & 
LGAS/SURFACE  reactions 

LIQU  ID/FOAM 

♦  c- phase  reactions 

♦  thermal  density  changes 

♦  mixture  properties 

♦  thermal  conduction 

♦  convection 

♦  molecular  diffusion 

♦  bubble  formation 


*{  MELTING 


UNREACTED  SOLID 

♦  thermal  conduction 

♦  convection 

♦  thermal  density  changes 


Figure  1.  Photograph  of  an  (a)  RDX-composite  propellant  (M43)  deflagrating  at  1.6  MPa  and  (b)  a  caricature  of  the 
3-phase  molecular  processes  involved. 
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Getting  now  to  the  formalism,  I  consider  a  semi-infinite  solid  combusting  in  the  steady  state  at 
constant  pressure.  The  spatial  coordinate,  taken  as  x,  extends  from  -oo  deep  in  the  solid  at  an 
initial  temperature  of  To,  to  the  solid/liquid  interface  at  -  xuq,  to  the  liquid/gas  surface  at  x  =  0,  and 
finally  to  the  region  of  equilibrium  gas  products  at  the  adiabatic  flame  temperature  Tf  at  x  =  oo  . 
Equations  conserving  mass,  atomic  species,  and  energy  must  be  solved  in  each  phase  subject  to  the 
boundary  conditions  for  that  phase.  The  intraphase  solutions  must  also  satisfy  the  equations  of 
continuity  at  each  phase  boundary  and  whatever  additional  constraints  are  imposed  by  the  surface- 
regression  mechanism. 

3.1  Intraphase  Conservation  Equations 

The  conservation  equations  within  each  phase  are  developed  and  discussed  elsewhere  (3),  but  will 
be  summarized  for  convenience  and  completeness  here.  The  +  and  -  superscripts  on  the  location 
superscripts  indicate  the  side  of  the  boundary  where  the  values  are  taken,  e.g.,  -x+Uq  means 

evaluated  at  the  liquid  side  of  the  interface  at  the  coordinate  x  =  -xLiq  . 


3.1.1  Solid  Phase 

Species  conservation: 

m^  =  d)kWk  k  =  1,2,..., iV  ,  (1) 

ax 


assuming  molecular  diffusion  is  negligible  on  a  combustion  timescale,  where  m  is  the  total  mass 
flux  and  an  eigenvalue  for  the  problem,  Yk  is  the  mass  fraction  of  the  kth  species,  d>k  is  the  net 
rate  of  production  of  species  k  due  to  chemical  reactions,  Wk  is  the  molecular  weight  of  the  kth 

species,  and  N  is  the  total  number  of  distinct  species  in  all  three  phases.  These  equations  are 
subject  to  the  domain  boundary  conditions 

Yk=Ykx  k  =  at  *  =  -».  (2) 

The  set  of  mass  fractions  {rp0}  reflect  the  composition  of  the  unreacted  propellant. 


Energy  conservation: 


d_ 

dx 


X 


Sol 


dT_ ' 

dX  y 


-  rite 


Sol 


dT_ 

dx 


NSo, 

£<w,  =o. 


(3) 


Both  the  thermal  conductivity  XSol  and  the  average  specific  heat  of  the  solid  mixture  cSol  are,  in 
general,  functions  of  the  independent  variable  x.  Wk  is  the  molecular  weight  and  hk  is  the 
enthalpy  of  species  k.  This  equation  is  subject  to  the  domain  boundary  conditions: 

T  =  T0  at  x  =  -oo  (4) 


and 
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T  =  Tm  at  x  =  -xLiq,  (5) 

where  T0  is  the  initial  temperature  of  the  propellant,  Tm  is  the  melting  point  of  the  propellant 
compound,  and  -xuq  is  the  coordinate  of  the  solid/liquid  boundary,  which  will  be  an  eigenvalue  of 
the  complete  combustion  problem.  It  should  be  noted  here  that  these  conservation  equations  for 
the  solid  phase  and  liquid  phase  in  the  next  section  embody  precepts  deriving  from  long  experience 
with  low-pressure  gas-phase  processes  and  may  not  be  as  general  as  presumed.  In  the  gas  phase  at 
sufficiently  low  densities  each  reaction  event  takes  place  in  essential  isolation.  Thermal  reaction 
coefficients,  obtained  by  averaging  over  velocities  and  reaction  cross  sections,  can  therefore  be 
used  to  characterize  reaction  events  anywhere  and  everywhere.  In  the  condensed  phase,  on  the 
other  hand,  reactive  events  do  not  take  place  in  isolation  but  in  a  dielectric  field  produced  by  the 
close  proximity  of  many  “spectator”  molecules.  This  dielectric  field  will  evolve  as  a  function  of 
the  instantaneous  configuration  of  all  the  molecules  in  the  vicinity  of  the  reacting  molecules  during 
the  course  of  the  reaction.  The  potential  therefore  exists  that  there  will  be  too  many  variables  to 
make  the  notion  of  thermal  rate  coefficients  an  adequate  and  useful  idealization.  Of  course,  at 
present  the  issue  is  moot  in  the  case  of  energetic  materials  as  the  reactions  are  virtually  unknown. 
An  MD  approach,  in  principle,  would  handle  these  complexities  in  a  natural  way. 

3.1.2  Liquid  Phase 

All  models  to  date  have  neglected  molecular  diffusion  in  the  liquid  phase  as  well,  but  I  retain  it  in 
the  general  formulation  here  because  I  later  show  that  it  may  be  of  some  importance. 

Species  conservation  in  liquid  phase: 

^-(pu,y,Vt)+m^-AtW^  0  k  =  1,2,...,  iV  ,  (6) 

UA  UA 

where  pLiq  is  the  liquid  density  and  Vk  is  the  diffusion  velocity  of  species  k.  These  equations  are 
subject  to  the  boundary  conditions: 

mYkXji*  +  k  =  1,2,...,N  at  x  =  -xLiq  (7) 

and 

Yk=Y[  k  =  l,2,...,N  at  x  =  0.  (8) 


Energy  conservation  in  liquid  phase: 


d_ 

dx 


( 


dT  \  . _u  dT 

-  -  me  - 

dx )  dx 


k 


2>,wa  =o. 


k= 1 


subject  to  the  boundary  conditions: 

T  =  Tm  at  x  =  -xLiq 


(9) 

(10) 
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and 


T  =  TS  at  x  =  0  . 


(11) 


3.1.3  Gas  Phase 

Species  conservation  in  gas  phase: 

d  dY 

-r(p0»W)  +  m-i-ffllW,=  0  *  =  1,2 . N.  (12) 

clx  ax 

subject  to  the  boundary  conditions: 

+  pLiqY[v[  =  mYf  +  pGXvf  k  =  1,2,...,N  at  x  =  0  (13) 


and 


dx 


=  0  k  =  l,2,...,N  at  x  =  co. 


Energy  conservation  in  gas  phase: 


d_ 

dx 


2  CJL 

^ Gas  , 

V  dx 


dT 


N 


N 


-  ^r-  Z  PoJXP,  -  Z  =  o, 

dx  V  ti 


subject  to  the  boundary  conditions: 


and 


T  =  T  at  x  =  0 


dT 

- =  0  at  x  =  oo. 

dx 


There  is,  finally,  the  mass  conservation  equation 


dm  _ 
dx 


(14) 

(15) 

(16) 

(17) 

(18) 


which  has  the  trivial  solution 

m  =  cons  tan  t  =  p''r,  (19) 

through  all  three  phases.  In  particular,  at  T0  ,  the  solid-phase  density  is  p°  and  the  linear  burning 
rate  is  r. 

3.2  Phase-Matching  Continuity  Conditions 

Some  of  the  above  boundary  conditions  are  expressed  in  quantities  that  are  unknown  at  the  outset 
and  coupled  to  solutions  in  adjacent  domains.  These  initially  unknown  quantities,  m  ,  Ts ,  xUq  , 
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and  {Yk°  ]  are  the  eigenvalues  for  the  complete  problem.  The  final  solutions  for  the  temperature 

and  mass  fractions  in  each  domain  must  satisfy  the  conservation  equations  with  their  respective 
boundary  conditions  in  each  phase  for  unique  values  of  the  eigenvalues.  In  order  to  accomplish 
this,  one  needs  further  constraints  on  the  problem.  These  are  to  be  found  in  the  inter-phase 
equations  of  continuity,  discussed  in  this  section,  and  in  the  surface-regression  mechanism, 
discussed  in  the  next  section.  Figure  2  shows  the  contributions  to  the  energy  fluxes  across  the 
phase  boundaries.  The  meaning  of  the  subscripts  and  superscripts  will  be  clear  from  the  figure 
contexts.  The  continuity  conditions  at  each  phase  interface  are  constructed  by  equating  the 
summed  contributions  on  each  side  of  a  particular  interface.  In  the  figure,  I  have  neglected 
contributions  from  the  kinetic  energy,  as  these  are  very  small  for  typical  propellant  burning  rates 
(of  order  0.01%  of  the  starting  enthalpies),  and  also  contributions  from  molecular  diffusion  in  the 
solid  phase,  because  it  will  be  too  slow  on  the  time  scale  of  importance  to  combustion.  The 
boundary  conditions  are  shown  in  Figure  2. 


Figure  2.  Energy  fluxes  at  the  phase  boundaries  arising  from  convection,  molecular  diffusion,  and  thermal 
conduction. 


3.2.1  Species-Flux  Continuity  at  Solid/Liquid  Boundary,  x  =  -xLiq 

This  relation  is  the  same  as  equation  7  and  thus  is  not  a  new  constraint. 
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3.2.2  Energy-Flux  Continuity  at  Solid/Liquid  Boundary,  x  =  -jc 


Liq 


-A, 


'Sol 


f  at\  Xliq  J*  f  at\ 

+  =-X 


\  dx  j 
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IT  \  -IM  N  N 
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3.2.3  Species-Flux  Continuity  at  Liquid/Gas  Boundary,  x  =  0 

This  relation  is  the  same  as  equation  13  and  thus  is  not  a  new  constraint. 

3.2.4  Energy-Flux  Continuity  at  Liquid/Gas  Boundary,  x  =  0 


“4 Liq 


CIL 

V  dx  j 


+  mYJYt°K'+  pu,Y,y°V°K  =  -Xam 


f  dT^ 

V  dx  j 


(20) 


+  m 


2X‘<*  +Pc„t,r"\Fhf  ,(2i) 


3.3  Surface  Regression  Mechanism 

In  general,  one  might  consider  that  many  mechanisms  contribute  to  surface  regression  during 
combustion:  evaporation  or  desorption  of  surface  species  without  change  of  identity,  reaction  of 
surface  species  to  produce  different  gas-phase  species,  reaction  of  gas-phase  species  with  surface 
species  to  produce  other  gas-phase  species,  and,  possibly,  physical  ejection  of  molecular 
aggregates  due  to  explosive  reaction  in  the  surface.  It  is  also  possible  that  some  combination  of 
these  mechanisms  might  be  active  in  the  condensed  phase  below  the  surface,  creating  bubbles 
which  are  subsequently  convected  to  the  surface  and  released.  Obviously,  it  will  be  a  very  great 
challenge  to  include  all  of  these  processes  in  a  full  combustion  model,  although  attempts  have  been 
made  to  model  RDX  using  a  mix  of  evaporation  and  bubble  formation  ( 4 ,  5).  In  these  treatments, 
one-dimensionality  was  preserved  through  the  artifice  of  continuous  porosity,  though  with 
unknown  accuracy.  In  my  opinion,  the  only  surface-regression  mechanism  that  has  been  treated 
with  reasonable  rigor  is  evaporation,  and  even  this  has  been  approximate  with  unknown  accuracy. 
The  remainder  of  this  section  will  be  devoted  to  a  discussion  of  the  approach  taken  to  an 
evaporative  mechanism  and  its  potential  shortcomings. 

3.3.1  Single- Component  Evaporation  Mechanism 

To  date,  the  only  evaporation  mechanism  used  in  published  burning-rate  model  developments  is 
based  on  the  following  reasoning.  Consider  a  pure  liquid  substance  in  equilibrium  with  its  vapor  at 
some  temperature  Ts.  I  will  refer  to  the  multiphase  presence  of  a  single  chemical  species  as  “single 
component.”  At  equilibrium,  the  mass  flux  impinging  on  the  surface  is  given  by  the  kinetic-theory 

result  ^nvW ,  where  n  is  the  molar  number  density  in  the  vapor,  v  is  the  mean  molecular  velocity 

at  Ts,  and  W  is  the  molecular  weight  of  the  species  being  considered.  If  a  is  the  fraction  of 
impinging  molecules  that  are  absorbed  into  the  surface  (i.e.,  the  accommodation  coefficient),  then 

the  mass  flux  actually  absorbed  into  the  surface  is  ^ anvW  .  Under  equilibrium  conditions,  the 
mass  flux  escaping  the  surface  equals  the  mass  flux  being  absorbed.  Thus,  the  escaping  mass  flux 
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can  be  expressed  in  terms  of  the  number  density  at  equilibrium  or,  through  the  equation  of  state, 
the  equilibrium  vapor  pressure  p.  Under  combustion  conditions,  the  vapor  pressure  at  the  surface 
is  less  than  the  equilibrium  value  due  to  depletion  by  reactions  and  molecular  diffusion.  Thus,  one 
can  express  the  net  regression  rate  of  the  surface  during  combustion  in  terms  of  the  equilibrium 
vapor  pressure  at  Ts,  the  total  pressure,  and  the  mole  fraction  of  the  mother-liquor  species  on  the 
gas-phase  side  of  the  surface,  X  .  These  arguments  are  summarized  graphically  in  Figure  3  and 
in  equation  22,  which  expresses  the  net  flux  of  evaporating  molecules  as  the  difference  between 
the  gross  escaping  flux  and  the  gross  condensing  flux,  i.e., 


m(xfr  ,T)  =  a 


W 


k2kRTsJ 


\TS)-Xl 


P total 
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Surface  Regression:  m  =  mescaping-  mcondeming 


Figure  3.  Single-component  evaporative  surface-regression  mechanism. 

3.3.2  Multicomponent  Evaporation  Mechanism 

On  the  face  of  it,  there  is  no  problem  with  the  rigor  of  the  arguments  just  presented.  The  difficulty 
comes  when  applying  it  to  even  the  simplest  combustion  problem.  Take  the  self-deflagration  of 
pure  frozen  ozone,  for  instance.  This  case  will  be  developed  more  fully  later,  but  it  serves  to 
illustrate  the  problem  with  single-component  evaporation.  Under  combustion  conditions,  by 
definition,  there  will  be  chemical  species  other  than  ozone  in  the  gas  phase.  These  will  be  products 
and  intermediates  of  ozone  decomposition,  i.e.,  CF  and  O.  These  other  species  will  also  impinge 
upon  and  be  absorbed  to  some  extent  in  the  surface.  Thus,  the  surface  will  have  a  multicomponent 
nature,  and  the  equilibrium  vapor  pressure  of  pure  ozone  will  no  longer  determine  the  mass  flux  of 
ozone  escaping  the  surface  but  also  be  a  function  of  the  mole  fractions  and  molecular  properties  of 
the  other  molecules  present.  In  addition,  the  mass  flux  leaving  the  surface  will  include 
contributions  from  the  surface-absorbed  products  and  intermediates  originating  from  the  gas-phase 
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reactions.  In  principle,  then,  the  single-component  evaporation  mechanism  is  intrinsically 
inappropriate  for  a  problem  involving  combustion.  Of  course,  as  an  approximation,  it  may  turn  out 
to  be  useful,  but  its  accuracy  is  impossible  to  assess  a  priori. 

3.4  Mathematical  Closure  of  the  3-Phase  Problem 

To  illustrate  the  final  posing  of  the  3-phase  mathematical  problem,  I  adopt  single-component 
evaporation  as  the  surface-regression  mechanism  and  assume  that  molecular  diffusion  in  the  liquid 
phase  can  be  neglected.  Neglect  of  molecular  diffusion  means  that  only  one  boundary  condition  in 
the  liquid-phase  domain  is  required  and  I  choose  this  to  be 

Yk=Yk~«*  k  =  1,2,...,N  at  x  =  -xLiq.  (23) 

—x~  —x+ 

Note  that  in  the  absence  of  diffusion,  Yk'Uq  =  Yk  L“‘ .  This  single  boundary  condition  replaces 

equations  7  and  8.  Neglect  of  diffusion,  through  replacement  of  equation  8,  also  greatly  simplifies 
the  problem  by  reducing  the  number  of  eigenvalues  to  three,  xUq  ,  Ts ,  and  m  .  One  begins  by 

providing  starting  estimates  for  these  three  eigenvalues  of  the  coupled  boundary- value  problems; 
call  them  xLiq ,  Ts ,  and  m  .  The  solid  phase  conservation  equations  arc  then  solved  using  the 

estimated  values  xLiq  and  m  .  This  solution  will  provide  values  of  j  Yk  Xuq  j  needed  by  the  liquid- 

phase  boundary  condition  equation  23.  After  solving  the  liquid-phase  conservation  equations  test 
to  see  if 

A4^F\  •  +'i‘Z!'r“'k",(7’J-Af‘(7’„r)].  (24) 

\  ax  I  ax )  i 

v  y  *=-%,■,  v  7  x=-xUq  * 

This  is  equation  20  rewritten  for  the  neglect  of  liquid-phase  diffusion.  The  quantity  in  brackets  is 
the  latent  heat  of  fusion  for  simple  substances.  If  the  equality  is  not  met  at  some  chosen  level  of 
approximation,  then  choose  another  value  of  xUq  and  repeat  until  equation  24  is  satisfied.  When 

equation  24  is  satisfied,  solve  the  gas-phase  conservation  equations  using  the  liquid-phase 
solutions  just  obtained  for  {  Y°  ],  then  test  to  see  if  the  surface  regression  mechanism  is  satisfied, 
i.e.,  if 

i 

4x°‘’f.)  ?  ?»«<]•  (25) 

-  2,tAT  J 

Normally,  equation  25  will  not  be  satisfied  at  this  point.  Choose  another  value  of  m  and  solve  the 
solid,  liquid,  and  gas  phases  again,  iterating  until  both  equations  24  and  25  are  satisfied.  Then  test 
to  see  if 
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This  is  equation  21  rewritten  for  the  neglect  of  liquid-phase  diffusion.  Normally,  equation  26  will 
not  be  satisfied  at  this  point,  so  choose  another  value  of  Ts  and  iterate  all  of  the  above  until 

equations  24-26  are  all  satisfied  to  the  required  degree  of  accuracy.  At  this  point,  one  will  have 
determined  the  eigenvalues  of  the  problem,  xuq,  Ts  ,  and  m  .  The  whole  process  can  be  more  easily 

visualized  as  a  logic  flow  chart  in  Figure  4.  Actually,  the  process  is  straightforward  and  successive 
guesses  for  the  eigenvalues  can  be  based  on  physical  reasoning  depending  on  the  results  of  each 
energy-flux  equality  test  so  that  convergence  of  each  eigenvalue  can  be  approached  with 
monotonic  decreasing  error.  The  solution  process  is  obviously  more  complicated  if  molecular 
diffusion  in  the  liquid  phase  is  not  neglected. 


4.  Models 


The  quantitative  modeling  of  the  burning  rate  of  solid  energetic  materials  really  began  in  the  dark 
hours  of  World  War  II  with,  predominantly,  the  efforts  of  Parr  and  Crawford  (6)  at  the  University 
of  Minnesota  and  Rice  and  Ginell  (7)  at  the  University  of  North  Carolina.  These  and  related 
wartime  works  were  published  in  a  special  issue  of  the  (then)  Journal  of  Physical  and  Colloid 
Chemistry.  At  that  time,  virtually  all  propellants,  gun  and  rocket,  used  some  combination  of  NC 
and  NG  as  their  energetic  ingredient.  Since  that  time,  composite  propellants  consisting  of 
ammonium  perchlorate  in  rubber  binders  (e.g.,  in  the  space  shuttle  boosters)  have  dominated  the 
solid  rocket  propellant  applications  with  solid  composites  based  on  HMX  in  polymeric  binders 
also  in  use.  Fielded  gun  propellants  are  still  dominated  by  nitrate-ester  propellants  although 
composites  of  these  conventional  propellants  with  NQ  are  common.  The  variety  of  propellant 
ingredients  exacerbates  the  difficulties  faced  by  combustion  modelers  both  because  of  their 
different  mechanisms  and  because  of  the  paucity  of  detailed  experimental  data  available  for  many 
of  them.  Recent  trends  are  toward  even  more  rapid  proliferation  of  new  chemical  ingredients  such 
as  oxetanes,  with  their  functional-group  tailorability,  and  azides,  with  their  attractive 
environmental  advantages.  In  addition  to  enhanced  performance  and  safety,  the  new  constraints  of 
minimal  environmental  impact  in  manufacture,  use,  and  demilitarization  are  now  driving  concerns 
in  the  development  of  new  propellants.  Coincident  with  the  emergence  of  many  promising  new 
energetic  materials,  with  attendant  dilution  of  experimental  characterization,  is  the  growing 
urgency  for  theoretical  guidance  in  the  formulation  of  propellants  incorporating  these  materials. 

It  is  generally  true  that  the  higher  the  performance  required  of  the  weapon  system,  the  smaller  the 
margin  of  safety  will  be  in  the  functioning  of  all  the  components  of  the  system  including  the 
propellant.  If  having  a  buming-rate  model  was  deemed  important  50  years  ago,  it  is  considerably 
more  so  with  today’s  new  mix  of  developmental  constraints  and  advanced  ingredients. 

In  this  section,  my  aim  is  not  to  give  a  comprehensive  history  of  buming-rate  modeling  but,  rather, 
to  provide  a  sense  of  the  conceptual  development  of  modeling  approaches.  This  background  is 
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Figure  4.  Flow  chart  illustrating  the  logic  for  determining  the  eigenvalues 
for  the  3-phase  problem  for  single-component  evaporation. 

essential  to  an  assessment  of  future  avenues  of  progress  in  the  field.  Three  subtopics  will  be 
addressed:  frozen  ozone,  RDX,  and  multi-ingredient  propellants.  This  progression  of  increasing 
system  complexity  allows  us  to  illustrate  some  of  the  detailed  mechanistic  challenges  facing  the 
model  builder  and  some  new  approaches  to  realizing  a  workable  tool  for  the  propellant  formulator. 
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4.1  Frozen  Ozone 


Frozen  ozone  is  the  simplest  chemical  system  falling  within  the  scope  of  3-phase  self- sustained 
deflagration.  Though  simple  from  a  theoretical  point  of  view,  it  is  anything  but  a  straightforward 
subject  for  experimental  investigation.  Its  propensity  to  detonate  is  legendary  and  the  attendant 
dangers  have  undoubtedly  inhibited  the  kind  of  extensive  measurements  of  burning  rate  that  one 
would  like  for  comparison  with  model  outputs.  On  the  other  hand,  a  wealth  of  high-quality 
experimental  data  has  been  obtained  on  thermophysical  properties  such  as  specific  heat,  thermal 
conductivity,  melting  and  boiling  points,  latent  heats,  reaction  paths  and  rates,  and  equations  of 
state.  This  comprehensiveness  and  reliability  of  the  input  database  on  frozen  ozone,  coupled  with 
its  simplicity,  makes  it  an  attractive  subject  for  modeling  despite  the  paucity  of  buming-rate  data. 
Its  conceptual  simplicity  encourages  and  enables  a  more  thorough  study  of  mechanisms  than  with 
any  other  energetic  material. 

Frozen  ozone  melts  at  about  80  K  and  has  a  normal  boiling  point  of  161.3  K.  The  rate  of  reaction 
in  the  condensed  phase  is  known  to  be  very  slow  compared  to  the  time  scale  of  self- sustained 
deflagration.  Thus,  all  of  the  uncertainties  of  describing  condensed-phase  reactions  are 
conveniently  (and  legitimately!)  sidestepped.  The  gas-phase  reaction  mechanism  is  known  with 
good  confidence  to  consist  of  the  following  three  reactions: 

O3  +  M  <->  O2  +  O  +  M  AH%8 15 k  =  +25.65  kcal  / mol ;  (I) 

03  +  O  <-»  02  +  02  AH°29815K  =  -93.41  kcal  /  mol ;  (II) 

O  +  O  +  M  <-»  02  +  M  A H°29815k  =  -119.06  kcal /mol .  (Ill) 


If  one  is  considering  a  pure  ozone  feedstock,  then  equation  22  describes  the  mass  flux  resulting 
from  evaporation.  The  only  measurement  of  the  linear  deflagration  rate  of  condensed-phase  ozone 
was  made  by  Streng  (8)  for  a  liquid  mixture  of  90%  O3  with  10%  CF.  The  total  mass  flux  leaving 
the  surface  in  this  multicomponent  case  is 

m  =  m(\-Y^  )+m0i,  (27) 


an  equation  first  used  by  Ben-Reuven  and  Caveny  (9)  in  connection  with  an  RDX  evaporation 
model.  m0  is  the  mass  flux  of  O3  alone.  Assuming  that  Y'(\  «  Y0'J  (i.e.,  no  liquid-phase  reactions 
or  molecular  diffusion),  this  expression  reduces  to 


m 


mo, 


(28) 


This  mass-fraction-equivalency  assumption  is  not  strictly  true  for  a  number  of  reasons  to  be 
discussed  subsequently,  but  it  may  be  an  adequate  approximation.  mQ  may  be  computed  using  a 

modification  of  equation  22  in  which  the  equilibrium  vapor  pressure  term  pe  is  replaced  by  the  O3 
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partial  pressure,  which  may  be  approximated  by  X0Of  p(0_  ,  an  expression  of  Raoult’s  law  for  ideal 
solutions.  X0Oi  is  the  O3  mole  fraction  at  the  surface,  which  is  assumed  to  be  approximated  by 
X0 “  ,  consistent  with  equation  28.  Of  course,  Raoult’s  law  is  only  an  approximation  that  can 
sometimes  result  in  considerable  error. 

An  example  of  a  calculation  (10)  of  the  burning  rate  of  pure  frozen  ozone  at  0.1  MPa  and  initial 
temperature  40  K  is  given  in  Figure  5.  As  indicated  in  the  figure,  the  computed  linear  burning  rate 
is  0.25  cm/s  at  0.1  MPa  and  an  initial  temperature  of  40  K.  By  comparison,  pure  RDX  burns 
~1  order  of  magnitude  more  slowly  at  a  sevenfold  higher  initial  temperature!  It  is  interesting  that 
the  surface  temperature  is  ~3  K  lower  than  the  boiling  point  at  this  pressure.  An  earlier  model  (11) 
assumed  that  the  surface  temperature  is  equal  to  the  boiling  point;  evidently,  in  this  case,  it  is  not  a 
bad  assumption.  Unlike  any  of  the  burning-rate  calculations  for  more  complex  systems,  such  as 
RDX,  the  calculations  by  Miller  (10)  include  the  effects  of  thermal  diffusion  in  the  gas  phase, 
although  the  burning  rate  is  changed  by  only  a  few  percent  if  this  process  is  neglected. 


Figure  5.  Computed  eigenvalues  and  profiles  for  the  steady-state  deflagration 
of  frozen  ozone  at  0.1  MPa  and  an  initial  temperature  of  40  K. 
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Another  mechanism  investigated  by  Miller  (10)  and  unique  to  this  ozone  study,  is  the  possibility  of 
a  heterogeneous  reaction  in  which  an  oxygen  atom  from  the  gas  phase  reacts  with  a  surface  ozone 
molecule  resulting  in  two  gas-phase  02  molecules: 

O(gas)  +  OsCliq)  — >  2  CFlgas)  AH°29S15K  =  -93.41  kcal/ mol .  (IV) 


This  reaction  apparently  has  never  been  measured  but  seems  not  only  plausible  but  probable.  It 
was  assumed  that  the  probability  is  unity  in  order  to  determine  the  maximum  effect.  Surprisingly, 
the  burning  rate  increased  only  1%  at  0.1  and  2.0  MPa,  despite  this  highly  exothermic  reaction 
occurring  at  the  surface,  where  it  adds  directly  to  the  heat  feedback.  In  addition  to  the  enhanced 
heat  feedback,  this  reaction  contributes  to  the  destruction  of  ozone  on  the  surface  and  should 
thereby  contribute  directly  to  the  regression  rate.  It  turns  out  that  the  O  atoms  near  the  surface  that 
are  consumed  by  the  heterogeneous  reaction  would  be  consumed  anyway  in  the  reaction  zone 
between  1  and  10  pm  from  the  surface  in  the  gas  phase,  and  reactions  this  close  to  the  surface 
contribute  their  heat  with  high  efficiency  to  the  heat  feedback  (Figure  6).  The  analysis  leading  to 
this  conclusion  is  based  on  an  important  concept  in  steady-state  combustion  that  can  be  gained 
from  a  study  of  the  previously  mentioned  conservation  equations.  Starting  with  equation  15  and 
assuming  that  the  mixture  thermal  conductivity  is  constant  and  that  the  specific  heats  for  all 
species  are  equal  and  constant,  one  can  show  (12)  that 


X 


Gas 


' dT ^ 
v  dx  )  X=Q- 


(29) 


where  q(x)  is  the  net  volumetric  rate  of  heat  release  in  the  gas  phase.  This  expression  indicates  that 

/t 

heat  released  within  a  characteristic  transport  distance,  Gf  ,  will  contribute  to  the  heat  feedback 

mcp 

with  high  efficiency.  Equation  29  is  valuable  in  very  complex  reaction  mechanisms  to  sort  out 
which  reactions  are  materially  affecting  the  burning  rate.  The  characteristic  distance  is  —10  pm  in 
Figures  5  and  6. 

The  ozone  problem  is  simple  enough  to  achieve  a  complete  and  unambiguous  analysis  of  the 
chemistry  in  each  part  of  the  flame.  The  monotonic  behavior  of  the  temperature  profile  in  Figure  5 
belies  the  underlying  complexity  of  parallel  and  sequential  reactions,  which  switch  directions 
during  the  course  of  the  flame,  as  seen  in  Figure  6.  A  curiosity  indicated  in  Figure  6  is  that 
spatially  the  first  reaction  step  in  the  ozone  decomposition  flame  is  the  production  of  ozone  by 
reaction  I  running  in  reverse.  In  more  complex  reaction  mechanisms  with  dozens  of  species  and 
hundreds  of  reactions,  it  is  impossible  to  completely  analyze  the  chemistry.  In  those  cases  sorting 
tools  such  PREAD  (13),  ChemPlot  (14),  and  EFEMAP  (15)  coupled  with  considerable  experience 
are  essential  to  gain  insights  into  the  mechanisms  affecting  both  the  burning  rate  and  the  flame 
features  such  as  the  dark  zone. 
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Source:  (10) 

Figure  6.  Rates  of  production  of  each  species  in  a  pure  frozen-ozone 
problem  source  at  0. 1  MPa  and  40  K  initial  temperature. 

The  f  and  b  suffixes  indicate  that  the  reaction  identified  is 
predominately  proceeding  in  the  forward  and  backward 
directions,  respectively.  The  whole  flame  is  divided  into 
zones  based  on  the  predominant  reactions  there. 

The  lone  measurement  of  a  condensed-phase  burning  rate  involving  ozone,  to  which  I  previously 
alluded,  was  for  a  liquid  mixture  of  90%  O3  and  10%  O2.  The  result  was  -0.4  cm/s;  our 
calculation,  utilizing  equation  28  and  Raoult’s  law,  produced  0.30  cm/s.  Sandri  (11)  has 
speculated  that  Streng’s  measured  rate  is  too  high  because  of  heating  of  the  liquid  by  the  flame  via 
conduction  through  the  containing-vessel  walls  and  by  radiation  from  the  flame.  Plausible  as  they 
might  be,  however,  the  degree  to  which  errors  from  these  sources  affected  the  burning-rate 
measurement  is  not  possible  to  determine.  Thus,  there  may  remain  a  significant  error  in  our 
buming-rate  model,  and  it  is  worthwhile  to  look  for  shortcomings  in  the  idealization  which  might 
account  for  the  discrepancy.  If  such  shortcomings  matter  to  the  ozone  case,  they  may  well  matter 
in  more  complex  systems  as  well. 

4.2  Deficiencies  in  the  Idealization 

4.2.1  Multicomponent  Evaporation 

As  pointed  out  in  section  3.3.2,  any  combustion  problem  driven  by  an  evaporative  surface- 
regression  mechanism  is  multicomponent  in  essence  since  product  species  will  coexist  in  the  liquid 
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surface  either  as  a  result  of  the  presence  of  more  than  one  ingredient,  condensed-phase  reaction  or 
molecular  diffusion  of  gas-phase  reaction  products  back  to  and  adsorption  onto  the  surface.  This 
issue  is  most  easily  discussed  within  the  context  of  the  O3/O2  mixture  problem.  Let  us  assume  that 
there  will  be  no  atomic  oxygen  in  the  liquid  surface  due  to  its  high  reactivity.  The  attractive 
intermolecular  forces  dominate  the  heats  of  vaporization,  and  the  strength  of  that  interaction  is 
ordered  as  follows:  O3-O3  >  O3-O2  >  O2-O2.  The  heat  of  vaporization  of  an  O3  molecule  from  a 
surface  of  both  O3  and  O2  molecules  will  therefore  be  less  than  that  from  a  surface  of  pure  O3 
because  the  latent  heat  is  just  a  measure  of  the  work  required  to  remove  an  O3  from  the  surface. 

On  the  other  hand,  the  heat  of  vaporization  of  an  O2  molecule  from  the  same  surface  mixture  will 
be  greater  than  that  for  O2  from  a  pure  O2  surface.  These  modifications,  in  turn,  alter  the 
equilibrium  partial  pressures  of  both  O3  and  O2  outside  the  mixture  surface,  and  the  equilibrium 
vapor  pressure  is  a  key  determinant  of  the  surface-regression  rate  by  equation  22.  Thus,  O3  will 
escape  at  a  faster  rate  and  O2  at  a  slower  rate  from  a  surface  mixture  of  the  two  substances  than 
would  be  expected  using  the  single-component  evaporation  formulation  coupled  with  equation  28. 


Beyond  the  greater  difficulty  of  computing  the  rates  of  escape  of  each  of  the  molecules  in  a  multi- 
component  description,  there  are  now  fundamental  consequences  to  the  posing  of  the  mathematical 
problem.  Though  the  mass  fluxes  of  O3  and  O2  leaving  the  surface  are  different,  in  order  to 
maintain  a  steady  state,  the  linear  rates  of  regression  of  both  species  must  be  equal,  i.e., 
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I  have  assumed  that  the  liquid  surface  is  composed  only  of  O3  and  O2.  This  means  that 


(34) 


iha  and  m0  can  be  obtained  from  equation  22,  modified  as  previously  described,  using  the  proper 

multicomponent  equilibrium  partial  pressures  for  each  species,  of  course.  Ideally,  one  would  also 
dispense  with  Raoult’s  law  by  computing  the  partial  pressures  based,  for  example,  on  model 
potential-energy  functions  as  described  in  section  5.1.  However,  I  am  now  left  with  a  new 
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unknown  or  eigenvalue  of  the  problem,  Y%  ,  with  equation  33  as  a  new  constraint.  I  now  must  add 

a  new  nested  loop  to  the  flow  chart  in  Figure  4.  While  this  is  a  serious  increase  in  computational 
complexity  for  the  ozone  problem,  it  becomes  an  overwhelming  increase  in  complexity  for  a 
compound  like  RDX  which  has  dozens  of  potential  species  in  the  surface.  Taking  the  general  case 
where  there  are  n  chemical  species  in  the  surface,  the  number  of  eigenvalues,  and  hence  nested 
loops,  will  be  ( n  +  2).  Solution  of  a  complex  set  of  equations  like  this  may  be  possible,  but  it  will 
undoubtedly  require  a  different  approach  than  nested  loops,  possibly  a  global  optimization  scheme 
in  which  trial  values  of  all  of  the  variables  are  changed  after  each  iteration. 

4.2.2  Liquid-Phase  Diffusion 

In  the  O3/O2  mixture  case,  the  inclusion  of  multicomponent  evaporation  will  mean  that  the  mass 
fractions  of  O3  and  O2  at  the  liquid  side  of  the  surface  are  different  from  the  proportions  in  the 
feedstock  at  -00.  This  is  illustrated  in  Figure  7.  Since  there  are  no  reactions  in  the  liquid  phase  (by 
justifiable  assumption)  the  discontinuity  in  feedstock/surface  mass  fractions  will  be  resolved  in 
nature  by  molecular  diffusion,  which  must  be  considered  as  an  attendant  consequence  of  treating 
multicomponent  evaporation.  Because  diffusion  is  generally  slow  in  liquids,  it  will  manifest  itself 
in  this  case  as  very  steep  gradients  in  the  species  just  below  the  surface.  Obviously,  a  whole  new 
class  of  supportive  data  will  be  required. 
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Figure  7.  Muticomponent  mixture  of  03  and  02  illustrating  how  the  differing 
rates  of  evaporation  of  the  two  components  result  in  depletion  of  02 
mole  fraction  and  enrichment  of  the  O3  mole  fraction  at  the  surface 
necessitating  the  consideration  of  liquid-phase  molecular  diffusion  to 
ensure  species  continuity. 
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4.2.3  Real-Gas  Equation  of  State 


Since  our  calculations  of  the  burning  rate  of  frozen  ozone  and  the  liquid  O3/O2  mixtures  assumed 
an  ideal-gas  equation  of  state,  it  is  worth  considering  if  real-gas  effects  play  a  role.  The  critical 
pressure  and  temperature  of  ozone  are  5.46  MPa  and  261.1  K,  so  the  shortfall  in  computed  burning 
rate  at  0.1  MPa  cannot  be  attributed  to  this  source.  However,  calculations  by  Miller  (10)  for  frozen 
ozone  were  performed  for  pressures  up  to  2  MPa.  At  2  MPa  and  40  K  initial  temperature,  the 
computed  surface  temperature  is  217  K.  Using  standard  tables  (16),  one  finds  that  ozone  gas  at  the 

pv 

surface  has  a  compressibility  factor  (^7  >  where  v  is  the  molar  volume)  of  about  0.75.  At  these 


conditions,  the  attractive  molecular  forces  are  dominating  and  number  densities  are  higher  than 
ideal,  speeding  up  the  reactions.  At  the  same  time,  for  pressures  near  the  critical  point,  the  heat  of 
vaporization  is  decreasing  rapidly;  this  means  that  a  given  amount  of  heat  feedback  will  vaporize 
more  surface  molecules.  Thus,  it  is  possible  that  use  of  a  real-gas  equation  of  state  could  have 
significant  impact  on  a  burning-rate  calculation.  To  our  knowledge,  there  are  no  published 
buming-rate  calculations  using  a  real-gas  equation  of  state. 


4.2.4  Phase  Separation 

Below  ~93  K,  mixtures  of  O3  and  O2  separate  into  two  phases,  the  upper  one  being  02  rich  and  the 
lower  one  O3  rich.  The  proportion  in  each  phase  depends  on  the  temperature.  At  the  boiling  point 
of  liquid  O2,  90.2  K,  for  example,  a  starting  mixture  of  50%  O3  by  weight,  O3  in  the  upper  layer 
ultimately  settles  to  12%  by  weight  and  38%  by  weight  in  the  lower  phase.  Streng’s  measurement 
of  the  burning  rate  of  a  90%  03/10%  O2  mixture  was  conducted  in  a  9-mm  Pyrex  tube  cooled  on 
the  outside  by  liquid  oxygen.  No  mention  is  made  of  the  presence  of  two  phases  during 
combustion  of  the  mixture,  but  if  the  temperature  were  actually  90.2  K,  two  liquid  phases  may 
have  been  present.  It  is  possible  that  this  phenomenon  could  explain  the  discrepancy  between  his 
measured  rate  and  our  computed  one,  since  the  02-rich  upper  layer  might  be  expected  to  evaporate 
faster.  On  the  other  hand,  this  would  lead  to  a  progressively  richer  mixture  as  the  O2  selectively 
escaped;  in  turn,  the  specimen  would  not  burn  at  a  steady  rate.  It  also  might  have  been  the  case 
that  the  mixture  surface  regressed  too  fast  to  establish  this  phase  separation.  Finally,  the  presence 
of  the  flame  inside  the  tube  may  have  warmed  the  liquid  mixture  a  few  Kelvins,  reaching  the  point 
where  O3  and  O2  are  miscible  in  all  proportions. 

4.3  RDX 

By  far,  the  most  intensive  efforts  to  compute  burning  rate  based  on  elementary  gas-phase  reactions 
have  been  directed  towards  neat  RDX.  It  is  known  that  RDX  chemically  decomposes  both  in  the 
solid  and  liquid  phases,  but  the  solid-phase  reaction  is  slow  and  probably  not  relevant  to 
combustion  time  scales  (17).  RDX  also  has  a  relatively  high  vapor  pressure  and  was  hypothesized 
to  gasify  during  steady  combustion  by  evaporation  first  by  Ben-Reuven  and  Caveny  (9),  whose 
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model  employed  global  reactions  in  the  gas  phase.  The  first  use  of  the  evaporation  mechanism  in  a 
model  with  elementary  gas-phase  reactions  was  by  Melius  (18)  ,who  also  assumed  the  single 
condensed-phase  reaction, 

RDX(Liq)  — »  3  CH20  +  3  N20.  (V) 

His  calculated  burning  rates  at  0.1  and  2  MPa  are  in  excellent  agreement  with  experimental  data. 

Other  elementary  gas-phase  reaction  RDX  models  combining  evaporation  and  liquid-phase 
reactions  include  those  of  Liau  and  Yang  (4),  Davidson  and  Beckstead  (5),  and  Prasad  et  al.  (19). 
The  last  of  these  works,  though  sharing  most  of  the  same  condensed-  and  gas-phase  reactions  as 
the  first  two,  is  a  semi-empirical  model,  requiring  the  experimental  value  of  the  surface 
temperature  to  obtain  the  burning  rate,  and  will  not  be  discussed  here.  Both  Liau  and  Yang  and 
Davidson  and  Beckstead  assume  that  reaction  V  occurs  and  one  other  liquid-phase  RDX 
decomposition  reaction.  In  the  Liau  and  Yang  model,  this  other  reaction  is 

RDX(Liq)  — »  3  HCN  +3/2  NO  +  3/2  N02  +  3/2  H20,  (VI) 

whereas  in  the  Davidson  and  Beckstead  model,  it  is 

RDX(Liq)  — »  3  H2CN  +3  N02.  (VII) 

Both  models  assume  that  RDX  can  form  bubbles  in  the  liquid  layer  and  that  the  following 
secondary  reaction  takes  place  in  the  bubbles: 

N02  +  CH20  — »  NO  +  CO  +  H20.  (VIII) 


The  reaction  parameters  chosen  by  each  of  the  authors  are  given  in  Table  1. 


Both  models  describe  the  melt  region  as  2-phase,  consisting  of  liquid  and  bubbles;  however,  since 
both  models  are  1-D,  this  2-phase  character  is  treated  in  the  conservation  equations  as  a  continuous 
“porosity,”  <f>,  defined  as  the  ratio  of  the  cross-sectional  area  occupied  by  the  gas  to  the  total  cross- 
sectional  area.  The  rate  of  evaporation  into  “bubbles,”  or  more  accurately,  into  the  porosity,  is 
computed  using  equation  22  at  the  local  subsurface  temperature  instead  of  the  surface  temperature 
and  the  “surface  area”  of  bubbles  defined  by  Liau  and  Yang  (4)  as 


AbHbUa=(36my,3t2/3  ,  (j)<]- 

Kbbles  =  (36m)1' 3 {l  -  (j>)2/3  ,  </>>y 


(35) 


where  n  is  the  bubble  number  density.  These  equations  were  given  without  explanation,  except  to 
say  that  n  was  to  be  determined  empirically.  Davidson  and  Beckstead  (5),  who  used  these  same 
equations,  evidently  understood  that  n  would  be  a  constant,  independent  of  pressure,  and  they 
chose  1  x  10  (13)  as  the  appropriate  value  without  further  explanation.  Liau  and  Yang  provide  no 
value  for  n  used  in  their  calculations. 
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Table  1.  Condensed-phase  reaction-rate  coefficient  parameters  assumed  by  different  models  for  RDX  in 
the  Arrhenius  form. 


Reaction 

V 

VI 

VII 

VIII 

A 

4.66  x  1018 

Melius  (18) 

B 

0 

E 

47.8  kcal/mol 

A 

6.0  x  1013 

1.6  x  1017 

802 

Liau  and  Y ang  (4) 

B 

0 

0 

2.77 

E 

36  kcal/mol 

45  kcal/mol 

13.73  kcal/mol 

A 

4.88  x  1011 

6.5  x  1014 

802 

Davidson  and  Beckstead  (5) 

B 

0 

0 

2.77 

E 

36  kcal/mol 

45  kcal/mol 

13.73  kcal/mol 

Notes: 

Burning-rate  calculations  based  on  these  models  are  given  in  Figure  8. 

The  Miller  (10)  3-phase  model,  assuming  no  condensed-phase  reactions,  is  applied  to  RDX,  and  its  results  are  also 

shown  in  Figure  8. 

A  is  in  appropriate  centimeter,  mole,  Kelvin,  and  seconds  units. 

The  consequences  of  these  varied  assumptions  and  data  on  the  computed  burning  rate  are  shown  in 
Figure  8.  For  comparison,  I  applied  our  3-phase  code  used  in  the  ozone  case  (10)  to  RDX  using 
the  same  input  data  and  reaction  mechanism  employed  by  Liau  and  Yang  with  the  exception  that 
no  condensed-phase  reactions  were  considered.  This  model  thus  assumes  that  surface  regression  is 
by  evaporation  at  the  surface  only,  driven  by  the  heat  feedback  from  reactions  in  the  gas  phase 
whereas  in  the  Liau  and  Yang  model  some  RDX  decomposes  in  the  condensed  phase  and  some 
RDX  evaporates  into  the  liquid-phase  porosity  (representing  bubbles).  The  close  agreement 
between  all  of  these  models  and  the  experimental  burning  rates  (20-26)  is  nothing  short  of 
astonishing.  Given  the  varied  inputs  and  assumptions,  one  is  tempted  to  conclude  that  the  single 
common  element  is  the  evaporative  mechanism  and  that  the  condensed-phase  reactions  are  not 
playing  a  significant  role  in  the  burning  rate.  However,  Davidson  and  Beckstead  use  a 
substantially  lower  vapor  pressure  (vapor  over  solid)  than  that  used  by  Liau  and  Yang  (vapor  over 
liquid)  and  find  a  low  sensitivity  of  the  burning  rate  to  vapor  pressure.  Moreover,  Davidson  and 
Beckstead  predict  a  fairly  large  amount  of  the  RDX  decomposes  in  the  liquid  state,  unlike  Liau 
and  Yang  (Table  2).  Unfortunately,  Melius  did  not  report  the  vapor-pressure  expression  used  in 
his  model,  and  his  results  are  not  otherwise  directly  comparable  to  the  other  models  because  his 
model  used  an  earlier  gas-phase  reaction  mechanism.  All  the  other  models  used  a  mechanism 
developed  by  Yetter  et  al.  (27)  based  on  Melius’s  (18)  but  extended  with  new  reactions  thought 
important. 
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Figure  8.  Several  3-phase  models  compared  to  experimental  RDX  burning 
rates. 

In  the  final  analysis,  no  definitive  comparison  is  possible  because  of  the  incomplete  reportage  of 
input  data  used  in  the  Melius  (18),  Liau  and  Yang  (4),  and  Davidson  and  Beckstead  (5)  models. 
Although  each  of  the  models  based  their  solution  of  the  gas-phase  conservation  equations  on  the 
PREMIX  code  (28)  developed  at  Sandia  National  Laboratories,  they  each  implemented  the  3- 
phase  numerical  solution  in  different  ways.  None  of  the  authors  described  these  methods  in 
sufficient  detail  to  enable  others  to  reconstruct  them  unambiguously.  One  is  left  with  the 
inescapable  conclusion  that  the  close  agreement  between  models  and  experiment  in  Figure  8  is  not 
as  astonishing  as  it  first  seems  but  the  result  of  finding  a  combination  of  input  parameters  that 
works.  In  the  end,  reconciliation  of  these  various  model  assumptions,  input  data,  and  numerical 
approaches  is  probably  not  productive  because  of  the  fundamental  and,  so  far,  irreducible 
uncertainties  associated  with  the  condensed-phase  processes,  principally,  but  not  limited  to,  the 
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Table  2.  Comparisons  among  RDX  models  of  selected  features  at  1  atm. 


Model 

Ts 

(K) 

Vapor  Pressure 
(atm)  at  600  K 

Sensitivity  of 
Burning  Rate  to 
Vapor  Pressure 

RDX  Fraction 
Decomposed  in 
Liquid  Phase 

(%) 

Melius  (18) 

549 

Not  reported 

Not  reported 

4.7 

Liau  and  Yang  (4) 

573 

0.69 

High 

<1 

Davidson  and  Beckstead  (5) 

595 

0.17 

Low 

25 

Miller  (10) 

633 

0.69 

High 

0 

chemical  reactions.  Modeling  the  condensed  phase  in  the  kind  of  detail  attempted  by  Liau  and 
Yang  and  Davidson  and  Beckstead  is  not  currently  supported  by  the  availability  of  adequate 
experimental  measurements  and  may  never  be. 

4.4  Multi-Ingredient  Propellant  Mixtures:  A  New  Approach 

In  view  of  the  situation  previously  described,  how  can  one  go  beyond  a  relatively  simple 
propellant  ingredient  like  RDX,  which  has  been  the  subject  of  the  most  extensive  modern  research 
to  date,  to  treat  multi-ingredient  propellants?  I  asked  this  question  of  myself  6  years  ago  and  in 
seeking  an  answer  had  an  idea  which  became  the  basis  for  launching  a  new  approach.  The  idea 
was  to  construct  a  hybrid-rigor  model  in  which  the  gas  phase  would  be  treated  in  full  elementary- 
reaction  detail  but  the  condensed  phase  and  surface  gasification  would  be  treated  in  semi-empirical 
fashion  using  the  Arrhenius-like  expression  first  used  by  Rice  and  Ginell  (7)  and  known  in  the 
propellant  community  as  the  “pyrolysis  law”: 

_  Es 

m(T )  =  As  e  RT'  ,  (36) 

coupled  with  a  set  of  hypothesized  decomposition  products  from  each  propellant  ingredient. 

The  decomposition  products  would  conform  to  proper  chemical  balance  and  be  selected  according 
to  available  experimental  results  and  theoretical  reasoning.  Furthermore,  each  propellant 
ingredient  would  be  assumed  to  decompose  without  interference  from  other  ingredients  (i.e., 
noninteractively).  A  full  development  of  this  model  for  nitrate-ester  propellants  is  given  elsewhere 
(29,  30),  but  the  basic  concepts  will  be  discussed  here. 

Of  course,  if  one  had  to  measure  a  pyrolysis  law  for  every  combination  of  ingredients,  there  would 
be  no  hope  of  predictive  capability  from  such  a  model.  However,  Anatoli  Zenin,  having  spent  the 
last  40+  years  refining  the  technique  of  measuring  surface  temperatures  of  numerous  propellant 
mixtures  and  single  propellant  ingredients  with  microthermocouples,  discovered  that  a  single, 
universal  pyrolysis  law  (23)  (Figure  9)  holds  for  a  wide  range  of  double -base  propellant 
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Gas-Phase  Species-Production  Rates 
in  Frozen-Ozone  Combustion 
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Figure  9.  Zenin’s  universal  pyrolysis  law  with  sampling  of  his  experimental  data 
(25)  for  a  wide  range  of  double-base  propellant  ingredient  proportions. 

formulations  with  NG  content  from  0%  to  50%  and  NC  with  percent  nitration  from  about 
11.5%  to  13.5%.  As  seen  in  Figure  9,  even  extreme  initial  temperature  data  are  reasonably  well 
accommodated  by  the  same  pyrolysis  law.  Furthermore,  the  same  double-base  pyrolysis  law  also 
works  for  additions  of  HMX  (23).  I  wondered  if  other  classes  of  ingredients  (e.g.,  nitramine- 
binder  systems)  might  display  this  same  universality.  If  so,  this  idea  might  lead  to  a  predictive 
tool,  at  least  for  members  of  these  classes.  To  test  this  hypothesis,  Zenin  (23)  applied  his 
embedded  microthermocouple  technique  to  a  wide  range  of  nitramine-oxidizer/polymeric-binder 
combinations  (Table  3).  Figure  10  shows  the  results  of  this  work  (31-33)  to  date.  First  of  all, 
notice  that  there  is  little  difference  between  RDX  and  HMX  materials  with  the  same  binder.  Also, 
most  binders  with  RDX  and  HMX  oxidizers  group  closely  about  the  pyrolysis  law  labeled 
RDXBA,  which  was  obtained  as  a  least-squares  fit  to  the  RDX/poly  3,3-bis(azidomethyl)  oxetane 
(BAMO)-poly  3-azidomethyl-3-methyl  oxetane  (AMMO)  data  alone.  It  is  of  great  interest  that  the 
CL20/PUNE  results,  while  closely  adhering  to  the  form  of  the  pyrolysis  law,  occupies  a  very 
different  region  on  the  graph.  HMX/PUNE,  on  the  other  hand,  is  in  excellent  accord  with  the 
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Table  3.  Zenin’s  nitramine/binder  test-material  formulations  (31-33). 


Test-Mixture  Designation 

Binder-Ingredient 

Proportions 

Oxidizer-Binder 

Proportions 

RDX  (or  HMX  )7B  AMO- AMMO 

50:50 

80:20 

RDX  (or  HMX)7BAMO-THF 

50:50 

80:20 

RDX  (or  HMX)/CBIH 

— 

80:20 

RDX  (or  HMX)/GAP1U 

— 

80:20 

RDX  (or  HMX)/GAP2 

— 

80:20 

CL20  (or  HMX)/PUNE 

20:80 

70:30 

a  50%  by  weight  of  particle  sizes  <  50  ps,  50%  in  range  150-300  pm,  99%  purity  of  RDX  and 
HMX. 

Notes: 

THF  =  tetrahydrofuran. 

CBIH  =  copolymer  of  butadiene  and  isoprene  with  hydroxyl  terminated  groups. 

GAP1U  =  glycidyl  azide-polyurethane  copolymer. 

GAP2  =  glycidyl  azide  polymer  (molecular  mass  of  2000). 

RDXBA  pyrolysis  law.  The  reason  for  this  difference  is  not  known;  perhaps,  the  rate-limiting  step 
in  the  decomposition  of  CL20  is  significantly  different  than  for  RDX  and  HMX.  The  tentative 
conclusion  from  these  measurements  is  that  propellants  composed  of  RDX  or  HMX  in  both  inert 
and  energetic  polymeric  binders  do  indeed  follow  a  universal  pyrolysis  law.  CL20  clearly  follows 
such  a  law  for  one  energetic  binder  and  may  for  other  binders  as  well.  Other  binders  are  currently 
being  investigated.  The  parameters  for  these  pyrolysis  laws  along  with  Zenin’s  double -base  law 
are  given  in  Table  4.  It  would  be  of  considerable  practical  interest  to  explore  the  theoretical  basis 
for  the  existence  of  these  universal  pyrolysis  laws;  they  seen  to  suggest  a  common  rate-limiting 
step  associated  with  the  condensed  phase  and/or  interfacial  region  for  each  family.  An 
understanding  of  this  behavior  could  well  lead  to  a  predictive  capability  for  pyrolysis  laws,  thereby 
removing  a  large  empirical  component  in  the  present  model. 

By  way  of  reinforcing  the  general  applicability  of  the  form  of  the  pyrolysis  law,  equation  36,  it  is 
worth  noting  that  neat  RDX,  HMX,  and  ammonium  dinitramide  (ADN)  all  may  be  well  described 
by  this  same  functional  relationship,  as  can  be  seen  in  Figure  11.  In  addition  to  the  experimental 
data,  I  have  plotted  in  Figure  1 1  the  results  of  a  least-squares  fitting  of  the  pyrolysis  law  to  results 
calculated  using  my  3-phase  code  for  both  frozen  ozone  and  neat  RDX.  It  can  be  seen  that  the 
relation  between  the  calculated  surface  temperatures  and  burning  rates  are  in  excellent  accord  with 
the  form  of  equation  36.  In  view  of  all  the  foregoing  evidence,  there  is  considerable  support  for 
the  general  applicability  of  the  pyrolysis  law,  and  this  should  encourage  efforts  to  calculate  the 
parameters  from  first  principles,  an  important  step  in  removing  the  empiricism  from  this  modeling 
approach  and  opening  the  door  to  much  wider  prediction  capability. 
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Figure  10.  Zenin’s  data  for  a  wide  variety  of  nitramine/binder  combinations  (31-33).  The  pyrolysis  law 
identified  as  NTRB  is  a  least-squares  fit  to  all  of  the  RDX/  and  HMX/binder  data.  The 
RDXBA  pyrolysis  law  is  a  fit  to  the  RDX/B AMO- AMMO  data  alone,  and  the  CL20PUNE 
pyrolysis  law  is  the  least-squares-fit  to  the  CL20/PUNE  data  alone. 
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Table  4.  Parameters  for  various  pyrolysis  laws,  equation  36. 


Designation 

As 

(g/cnf-s) 

Es 

(cal/mol) 

DB 

1800 

9935 

RDXBA 

10470 

1371 

NTRB 

2004 

11827 

CL20PUNE 

1.868  x  106 

16032 
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Figure  11.  Pyrolysis  laws  for  a  number  of  neat  energetic  materials.  The  RDX  Miller  calculation  was  performed 

with  the  3-phase  model  previously  used  for  frozen  ozone  using  the  same  input  data  for  RDX  as  Liau  and 
Yang  (with  the  exception  of  no  condensed-phase  reactions).  All  data  and  calculations  are  for  an  initial 
temperature  of  293  K  except  for  frozen  ozone,  which  was  at  40  K. 

As  previously  indicated,  the  model  requires  the  identity  and  mole  fractions  of  the  chemical  species 
resulting  from  the  condensed-phase  processes  (i.e.,  the  gases  leaving  the  surface).  These  species 
will  be  far  from  the  equilibrium  distribution,  and,  therefore,  one  has  only  chemical  balance  to 
constrain  the  set  of  surface  products  absolutely.  With  ingredients  as  complex  as  those  for 
propellants,  this  means  that  there  are  usually  very  many  possible  product  sets  to  consider. 
Experimental  thermal  decomposition  studies  may  be  useful  guides  to  selecting  an  appropriate 
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product  set  but  not  perfect  since  those  experiments  may  not  accurately  mimic  the  conditions  of  the 
burning  surface  (temperature  and  heating  rate)  and  probably  all  suffer  from  significant  unknown 
changes  in  the  surface  products  as  a  result  of  very  fast  reactions  close  to  the  surface.  Use  of 
general  theoretical  reasoning  can  also  help  limit  the  number  of  probable  product  sets.  In  the  end, 
however,  an  important  selection  criterion  must  be  whether  the  proposed  product  set  results  in  a 
calculated  burning  rate  close  to  the  experimental  one.  Use  of  this  criterion,  of  course,  assumes  that 
the  pyrolysis  law,  the  gas-phase  reaction  mechanism,  and  all  the  other  input  data  the  model  needs 
are  accurate.  The  predictive  ability  of  the  model  for  the  whole  propellant  is  redeemed  by  insisting 
that  the  decomposition  products  from  each  ingredient,  “calibrated”  ideally  for  that  ingredient 
alone,  remain  consistent  for  all  propellant  formulations  using  that  ingredient.  This  strategy  is 
probably  the  best-compromise  approach,  but  one  should  recognize  that  the  assumption  of  non¬ 
interactivity  in  decomposition  might  be  a  better  approximation  for  some  systems  than  others. 

Since  all  fielded  gun  propellants  include  at  least  some  NC  and  certain  propellants,  such  as  M10, 
have  98%  NC,  any  burning-rate  model  of  practical  importance  will  have  to  deal  with  NC  as  an 
ingredient.  Not  only  is  NC  a  complex  long-chain  polymer,  but  the  repeat  units  vary  among  four 
different  types  depending  on  whether  the  glucose  ring  is  triply,  doubly,  singly  nitrated  or 
unnitrated  altogether.  Figure  12  shows  a  triply  nitrated  repeat  unit.  NC  is  characterized  for 
propellant  use  by  its  average  percent  nitrogen  (%N)  by  weight,  and  this  quantity  varies  typically 
from  -11.5  to  -13.5  %N.  Pure  cellulose  trinitrate,  cellulose  dinitrate,  and  cellulose  mononitrate 
have  14.14,  11.11,  and  6.76  %N,  respectively.  The  work  of  Leider  and  Seaton  (34)  showed  that  a 
12.3  %N  NC  has  all  three  nitration  states  present  but  no  unnitrated  sites.  We  developed  a  Monte 
Carlo  model  (29,  30)  to  determine  for  an  NC  specimen  of  given  percent  nitrogen  what  is  the 
distribution  among  these  three  nitration  states.  Predictions  of  this  model  are  compared  to  the 
nuclear  magnetic  resonance  (NMR)  data  of  Todd  and  Glasser  (35)  in  Figure  13.  The  burning-rate 
model  then  assumes  that  NC  consists  of  three  separate,  hypothetically  pure  ingredients 
corresponding  to  cellulose  tri-,  di-,  and  mononitrate  in  the  proportion  indicated  by  the  Monte  Carlo 
model.  The  densities  of  these  hypothetically  pure  nitrates  are  obtained  from  a  MD  calculation  (36) 
using  the  COMPASS  force  field. 

The  model  handling  of  multi-ingredient  propellants  is  summarized  graphically  for  JA2  in 
Figure  14.  This  figure  does  not  represent  the  current  decomposition-product  sets  but  serves  to 
illustrate  the  principles  behind  the  method.  The  CYCLOPS  (29,  30)  code,  named  for  the  mythical 
race  of  creatures  who  forged  thunderbolts  for  Zeus,  is  the  computer-code  embodiment  of  the 
buming-rate  model  under  discussion.  CYCLOPS  accepts  the  weight  percents  of  each  ingredient 
and  checks  to  see  if  NC  is  among  them.  If  so,  then  the  Monte-Carlo  subroutine  is  called  and  the 
nitrate  state  distribution  determined.  The  ingredient  list  is  then  expanded  to  include  the  three  NC 
subingredients  (cellulose  tri-,  di-,  and  mononitrate).  The  condensed-phase  decomposition 
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Figure  12.  Triply  nitrated  NC  repeat  unit. 

products  for  each  of  the  ingredients  on  the  expanded  list  are  then  obtained  from  the 
decomposition-products  database,  and  the  net  set  of  products  computed  according  to  the 
proportion  of  each  ingredient.  A  trial  value  for  the  linear  burning  rate  is  then  read  in  from  the 
problem  input  file  and  converted  to  mass  flux  using  the  code-computed  propellant  density.  The 
surface  temperature  for  this  mass  flux  is  then  computed  from  the  pyrolysis  law,  and  these  values 
are  passed  to  a  modified  version  of  the  PREMIX  (28)  code,  which  is  called  as  a  subroutine  to  solve 
the  gas-phase  conservation  equations.  From  this  solution  the  heat  feedback  to  the  surface  is 
computed  and  compared  with  the  heat  feedback  required  to  transform  the  propellant  ingredients  at 
their  initial  temperature  to  the  condensed-phase  decomposition  products  at  the  trial  surface 
temperature  using  the  trial  mass  flux.  New  trial  values  are  automatically  provided  by  the  code 
until  the  two  heat  fluxes  are  in  satisfactory  agreement.  The  criterion  for  convergence  is  actually 
expressed  in  terms  of  acceptable  changes  in  successive  guesses  of  the  mass  flux,  which  is  the 
unknown  sought. 
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Substitution  Level 


Figure  13.  Monte  Carlo  model  for  distribution  of  repeat  units  among  cellulose  tri-,  di-,  and 

mononitrates  for  an  NC  specimen  of  given  percent  nitrogen.  Comparison  is  made  of  the 
model  predictions  with  the  NMR  data  of  Todd  and  Glasser  (55). 
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Figure  14.  Conceptual  deconstruction  of  a  propellant  containing  NC  into  subingredients  and  then  into  net 
condensed-phase  decomposition  products  entering  the  gas  phase. 

Results  for  the  burning  rate  of  M10  and  JA2  are  shown  in  Figures  15  and  16  compared  to  the 
experimental  data  of  Miller  (22),  Juhasz  (37),  and  Atwood  et  al.  (38)  for  M10  and  Miller  (22)  and 
Juhasz  et  al.  (39)  for  JA2.  An  extensive  analysis  (29,  30)  of  the  chemical-kinetic  origin  of  the 
inflection  in  the  M10  curve  has  been  made.  It  was  discovered  that  those  reactions  controlling  the 
dark-zone  length  at  low  pressures  (<10  MPa)  have  little  influence  over  the  burning  rate  at  low 
pressures,  but,  to  our  surprise,  a  major  effect  on  the  burning  rate  at  high  pressures,  where  a  dark 
zone  does  not  even  exist.  Thus,  studies  of  the  dark-zone  chemistry  experimentally  accessible  at 
low  pressures  actually  are  probing  those  reactions  with  critical  influence  over  the  burning  rate  at 
high  pressures.  This  finding  provides  new  impetus  to  experimental  dark- zone  investigations, 
especially  for  nitramine  propellants  because  their  dark-zone  chemistry  is  less  well  known  than  that 
of  nitrate-ester  propellants.  How  well  CYCLOPS  predicts  the  structure  of  the  dark  zone  is  shown 
in  Figure  17,  where  the  temperature  and  NO  profiles  are  seen  to  compare  very  favorably  with 
experiment.  Table  5  gives  a  comparison  of  predicted  and  measured  major  species  for  a  double¬ 
base  propellant  similar  to  M9;  again,  the  agreement  is  excellent.  The  experimental  measurements 
used  for  comparison  are  those  of  Heller  and  Gordon  (40),  Lengelle  et  al.  (41),  and  Vanderhoff  et 
al.  (42),  and  they  were  renormalized  for  direct  comparability  by  Vanderhoff  et  al.  (42). 
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Figure  15.  Comparison  of  CYCLOPS-code  calculations  of  burning  rate  of  M10 
propellant  with  experimental  data. 

A  final  example  is  given  of  more  recent  work  (43)  on  RDX  in  an  energetic  thermoplastic  elastomer 
(ETPE)  binder.  The  excellent  comparison  with  experimental  buming-rate  data  shown  in  Figure  18 
is  tempered  by  the  poor  prediction  of  the  thermal  structure  in  Figure  19.  This  inconsistency  might 
at  first  be  thought  to  suggest  that  the  good  buming-rate  prediction  is  simply  fortuitous,  but  this  is 
not  necessarily  the  case.  As  we  learned  in  the  case  of  M10,  the  burning  rate  and  the  dark- zone 
length  at  low  pressure  may  be  controlled  by  different  sets  of  reactions,  so  it  is  possible  that  the 
good  burning-rate  prediction  and  poor  thermal  structure  prediction  may  be  due  to  imperfections  in 
the  chemical  mechanism  describing  the  nitramine  dark-zone  chemistry.  In  all  propellants  with 
RDX,  we  assume  that  the  RDX  evaporates  unchanged,  just  as  I  had  assumed  for  neat  RDX; 
however,  binder  ingredients  are  assumed  to  decompose  in  the  condensed  phase. 
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Figure  16.  Comparison  of  CYCLOPS-code  calculations  of  burning  rate  of  JA2 
propellant  with  experimental  data. 

The  foregoing  examples  give  evidence  of  the  promise  of  the  CYCLOPS  code  in  providing  both 
reasonable  predictions  of  burning  rate  as  well  as  details  of  the  gas-phase  flame  structure  and 
detailed  chemistry.  At  present,  the  code  is  limited  to  certain  families  of  propellant  ingredients,  but 
the  number  of  families  will  likely  increase  as  other  systems  are  studied.  However,  perhaps  more 
important  is  the  growing  probability  that  one  might  be  able  to  predict  pyrolysis  laws  from  first 
principles  considerations.  Certainly,  greater  theoretical  capabilities  to  predict  the  final 
condensed-phase  decomposition  products  of  propellant  ingredients  can  be  brought  to  bear.  Thus, 
the  CYCLOPS  code  may  well  provide  the  best  computational  vehicle  for  incorporating  results  of 
new  theoretical  approaches  to  the  condensed-phase  and  surface  processes.  How  this  might  be 
accomplished  while  maintaining  mathematical  tractability  of  the  code  is  suggested  in  section  5. 
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Figure  17.  Comparison  of  CYCLOPS-code  predictions  of  dark-zone  thermal 

structure  compared  with  the  experimental  data  of  Vanderhoff  et  al.  (42) 


Table  5.  Comparison  of  major  species  mole  fractions  in  the  dark  zone  of  double-base  propellant  with  various 
experimental  measurements. 


Parameters 

Heller  and 
Gordon3 

Lengelle  et  al.b 

Vanderhoff  et  al.c 

CYCLOPS 
(Present  Model) 

P  (MPa) 

1.6 

0.9 

1.7 

1.7 

Dark-zone  temperature  (K) 

1600 

1500 

1500 

1543 

NO 

0.24 

0.21 

0.24 

0.25 

CO 

0.33 

0.38 

— 

0.32 

h2 

0.08 

0.08 

— 

0.08 

n2 

0.04 

0.02 

— 

0.004 

h2o 

0.20 

0.20 

— 

0.19 

co2 

0.10 

0.09 

— 

0.10 

HCN 

0.004 

— 

— 

0.004 

CH4 

0.008 

0.026 

— 

0.009 

c2h4 

0.008 

0.008 

— 

0.001 

“  Source:  (40). 
b  Source:  (41). 
c  Source:  (42). 
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Figure  18.  Comparison  of  CYCLOPS-code  calculations  of  the  burning  rate  of  and 

RDX/thermoplastic  elastomer  (TPE)  propellant  using  two  different  pyrolysis 
laws  with  the  experimental  data  of  Zenin  (32). 


Figure  19.  Comparison  of  CYCLOPS-code  predictions  of  dark-zone  thermal  structure 
for  an  RDX/TPE  propellant  with  the  microthermocouple  data  of  Zenin  (32). 
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In  summary,  the  CYCLOPS  code  is  based  upon  a  few  simple  propositions: 

•  The  surface  regression  may  be  adequately  idealized  as  1-D.  This  condition  is  probably  met 
for  most  homogeneous  propellants.  It  may  also  be  approximately  true  for  many  composite 
propellants  where  the  burning  rate  is  not  a  strong  function  of  oxidizer  particle  size.  Surface 
melt  layers,  for  example,  may  provide  effective  premixing  of  ingredients  where  the  oxidizer 
particle  sizes  are  not  too  large. 

•  The  surface  regression  may  be  described  by  a  pyrolysis  law,  equation  36.  I  have  shown  that 
this  condition  has  been  met  by  a  large  number  of  neat  energetic  materials  and  mixtures  of 
typical  propellant  ingredients. 

•  The  overall  products  of  condensed-phase  decomposition  may  be  estimated  with  sufficient 
accuracy.  Chemical  balance  and  results  of  thermal  decomposition  experiments  are  sources 
of  guidance  here,  though  only  the  former  can  be  considered  as  an  absolute  constraint  because 
of  the  high  probability  of  secondary  reactions  in  decomposition  experiments. 

•  The  decomposition  of  one  propellant  ingredient  does  not  affect  the  decomposition  of  other 
ingredients,  i.e.,  the  decomposition  of  a  multi-ingredient  propellant  may  be  described  as  the 
superposition  of  the  independent  decompositions  of  each  of  its  ingredients.  This  may  well  be 
a  better  approximation  for  some  ingredients  than  for  others.  Short  of  describing  the  mixture 
decomposition  theoretically  on  a  molecular  scale,  there  is  little  recourse  to  this  assumption. 
However,  it  is  often  found  that  significant  changes  in  detailed  combustion  mechanisms  have 
a  remarkably  weak  influence  on  the  burning  rate,  that  quantity  being  a  highly  integrated 
consequence  of  myriad  underlying  details.  CYCLOPS  exploits  this  relative  insensitivity  to 
mechanistic  details. 

Finally,  since  the  CYCLOPS  code  is  the  first  code  of  its  kind  (though  still  undergoing 
development  to  improve  robustness  for  general  use),  one  can  anticipate  ways  in  which  it  might  be 
put  to  practical  use  even  at  its  current  stage  of  development.  First,  the  CYCLOPS  code  can  be 
used  to  optimize  ingredient  proportions  to  achieve  a  target  burning  rate.  Also,  CYCLOPS  could 
be  used  to  explore  subtle,  previously  unexplained  effects  of  formulation  on  burning  rate.  For 
example,  German-made  JA2  bums  about  20%  faster  than  U.S.-made  JA2  with  identical  ingredient 
proportions .  However,  the  two  differ  by  different  specifications  on  the  NC  component.  While  the 
average  percent  nitrogen  for  the  two  propellants  is  the  same,  the  two  materials  are  made  from 
blends  of  two  lots  of  NC  with  different  percent  nitrogen.  Because  our  nitrate- state-distribution 
Monte  Carlo  code  predicts  that  the  surface  products  for  these  two  propellants  will  be  different,  it 
will  be  worthwhile  to  see  if  CYCLOPS  can  predict  the  burning-rate  difference.  Further, 
CYCLOPS  could  be  used  to  help  set  manufacturing  tolerances  in  the  specifications  of  ingredient 
proportions  in  military  propellants.  The  CYCLOPS  code  could  be  used  to  determine  how  any 
given  set  of  tolerances  will  map  into  variations  of  the  burning  rate  and  these  burning-rate 
variations  could  be  judged  by  using  them  as  inputs  to  interior-ballistic  codes.  The  burning  rate  at 
low  pressure  has  been  shown  to  be  sensitive  to  the  heats  of  formation  of  propellant  ingredients  in 
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certain  cases  (29,  30).  To  the  extent  that  ingredient  purity  is  reflected  in  the  heat  of  formation 
(e.g.,  due  to  dinitroglycerine  impurity  in  trinitroglycerine),  the  CYCLOPS  code  could  be  used  to 
determine  how  ingredient  purity  affects  the  burning  rate,  thereby  enabling  rationalization  of 
quantitative  tolerances  for  these  ingredients.  Finally,  the  CYCLOPS  code  could  determine  the 
effect  of  chemical  modifiers  on  the  burning  rate  without  mixing  chemicals.  Though  not  replacing 
the  need  to  mix  and  test,  using  the  code  as  a  screening  agent  and  test  bed  for  ideas  could  improve 
efficiency  and  possibly  suggest  new  compounds  to  test. 


5.  Challenges  and  Opportunities 


The  semi-empirical  aspects  of  the  CYCLOPS  code  need  be  tolerated  only  until  theoretical 
advances  obviate  their  necessity.  The  next  level  of  improvement  may  well  come  from  MD 
descriptions  of  the  condensed  phase  and  surface  phenomena.  However,  this  more  fundamental 
level  of  treatment  will  have  its  own  set  of  approximations  and  limitations.  In  order  to  deal  with 
condensed-phase  reactions,  for  example,  considerable  progress  will  have  to  be  made  in  the 
parameterization  of,  and  experience  with,  reactive  force  fields.  In  my  judgment,  one  will  never 
want  to  treat  the  gas  phase,  with  its  dozens  of  species  and  hundreds  of  reactions,  in  this  way.  To 
do  so  would  discard  more  than  50  years  of  hard-won  kinetics  research  gains.  One  will  therefore  be 
faced  with  merging  the  discrete  and  continuous  descriptions  into  a  single,  tractable  mathematical 
entity.  I  believe  that  this  will  be  no  easy  task,  as  the  MD  approach  is  inherently  time-dependent 
with  (computationally  intensive)  explicit  statistical  averaging  and  the  steady- state  continuum- 
mechanics  approach  outlined  earlier  in  this  work  is  inherently  time-independent  with  implicit 
statistical  averaging.  The  most  fruitful  approach  may  be  to  develop  simple,  idealized  continuum 
submodels  calibrated  to  stand-alone  MD  models  of  a  limited  set  of  phenomena.  For  example,  a 
reactive-MD  model  of  the  condensed  phase  would  provide  a  set  of  condensed-phase 
decomposition  products  needed  by  the  CYCLOPS  code.  Similarly,  an  MD  model  of 
multicomponent  evaporation  and/or  pyrolysis  could  be  used  to  determine  a  pyrolysis  law  for 
CYCLOPS.  Results  of  the  continuum  submodels,  depending  on  their  scope,  could  then  be  easily 
incorporated  into  either  a  3-phase  model  or  a  code  such  as  CYCLOPS,  retaining  mathematical 
tractability  and  efficiency.  Since,  prior  to  the  development  of  CYCLOPS,  I  had  made  a  start  on 
incorporating  multicomponent  evaporation  into  our  3-phase  code,  it  will  be  used  as  an  example 
here  of  the  proposed  approach. 

5.1  “Molecular”  Continuum  Model  of  Multicomponent  Evaporation 

I  consider  that  a  molecule  evaporating  from  a  liquid  surface  (Figure  20)  is  at  some  distance  cl 
above  the  surface.  I  further  assume  that  the  molecule  follows  a  tubular  path  from  deep  within  the 
liquid  to  points  outside  the  liquid;  the  tube  has  radius  ao .  The  conception  here  is  that  the  radius  of 
the  tube  is  of  the  order  of  the  “radius”  of  the  escaping  molecule,  and  what  is  a  tortured  path 
geometry  in  reality  is  idealized  as  straightened  to  a  cylindrical  path.  Since  the  escaping  molecule 
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Figure  20.  Model  for  continuous-phase  molecular  forces  experienced  by  a  molecule  evaporating  from  a  liquid  surface. 

is  at  a  distance  r  from  all  the  molecules  in  the  differential  volume  element,  dV  =  2n  adadz ,  the 
potential  energy  of  interaction  is  dO  =  (f){r)n  dV.  where  n  is  the  number  density  of  liquid 
molecules  (Figure  20).  If  I  assume  a  Lennard-Jones  (L-J)  interaction  potential,  where  s  is  the 
well  depth  and  a  is  the  collision  diameter,  i.e., 
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and  integrate  over  all  the  molecules  in  the  liquid,  one  obtains 
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Two  important  extensions  can  be  made  to  this  formulation.  The  first  is  to  add  the  effects  of 
interaction  between  the  escaping  molecule  and  the  gas-phase  molecules.  This  modifies  the 
potential  to 
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The  other  important  extension  is  to  the  multicomponent  case.  This  requires  a  sum  over  the 
interactions  between  the  escaping  molecule  of  species  i  with  the  other  species  j,  of  total  different 
kinds,  N.  Thus,  the  multicomponent  version  of  the  interaction  potential  between  an  evaporating 
molecule  of  species  i  at  a  distance  cl  from  the  surface  and  all  other  species  both  in  the  liquid  and  in 
the  gas  phase  is 
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with  its  ancillary  definitions  in  equations  39-41.  These  equations  may  seem  awkward  but  are 
easily  and  quickly  evaluated  in  a  computer  code  and  represent,  in  a  very  general  way,  enormously 
complex  phenomena.  For  example,  they  describe  quantitatively  how  the  heat  of  vaporization  goes 
to  zero  at  the  critical  point,  where  the  liquid-  and  gas-phase  densities  become  indistinguishable. 


A  necessary  test  of  the  model  is  the  accuracy  with  which  the  previous  equations  predict  the  heats 
of  vaporization  of  pure  substances.  The  heat  of  vaporization  is  obtained  by  taking  the  limit  in  the 
previous  equations  as  cl  64  in  equation  38.  The  result  is 
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Application  of  this  model  is  made  to  61  data  sets  for  polar  and  nonpolar  molecules  using  the 
compilation  of  L-J  parameters  in  Reid  and  Sherwood  (44)  and  Reid  et  al.  (45).  Calculation  results 
shown  in  Figure  21  are  obtained  by  performing  a  non-linear  least-squares  fit  of  the  heat-of- 
evaporation  experimental  data  to  equation  44  using  Pe,  defined  as  follows,  as  an  adjustable 
parameter  to  mediate  between  the  raw  theoretical  result  Hcaic  and  the  best  estimate  value  Hest\ 

Hest=PeHcalc.  (45) 


Assuming  that  — 


is  unity,  the  best-fit  value  of  Pe  turns  out  to  be  1.10.  Its  proximity  to  unity 


suggests  that  the  physical  basis  of  model  assumptions  are  reasonable.  The  standard  deviation  of 
the  error  using  the  optimized  value  of  Pe  is  about  ±12%.  This  error  is  surprisingly  small 
considering  that  the  L-J  potential  is  not  generally  as  faithful  as  an  exponential-6  potential  for 
nonpolar  molecules  and  is  has  even  worse  fidelity  for  polar  molecules,  which  are  abundantly 
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Figure  21.  Accuracy  of  the  simple  heat-of-vaporization  theory  using  61  L-J 
parameter  sets  for  both  polar  and  nonpolar  molecules. 


represented  in  the  data  set.  In  fact,  even  without  the  adjustment  parameter  Pe,  the  standard 
deviation  of  the  prediction  of  Hvap  is  about  16%.  This  good  performance  might  be  improved  upon 
by  either  using  another,  better  interaction  potential  or  by  possibly  modifying  the  model  to 
recognize  the  discrete  nature  of  close-neighbor  molecular  interactions. 

By  itself,  the  previously  mentioned  model  is  not  sufficient  to  treat  the  multicomponent  evaporation 
problem.  One  needs  also  to  describe  the  dynamics  of  the  evaporative  escape  mechanism. 

Drawing  upon  kinetic  theory  to  construct  the  flux  of  molecules  leaving  the  liquid  surface,  rout , 

one  obtains 


(46) 


where  the  quantity  in  parenthesis  is  the  density  of  molecules  with  energies  greater  than  Hvap,  and 
ve  is  the  average  velocity  of  those  molecules  with  enough  energy  to  escape  the  heat-of- 

vaporization  barrier.  Assuming  that  the  molecules  in  the  surface  are  equilibrated  at  the  surface 
temperature,  Ts,  it  may  be  shown  that  the  average  escape  velocity  is 
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where  W  is  the  molecular  weight  of  the  escaping  molecules  and  the  minimum  velocity  for  escape 
Ve  IS 
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(49) 


The  equilibrium  vapor  pressure  is  determined  by  equating  the  outward  and  inward  fluxes  at 
equilibrium  (see  arguments  of  section  3.3.1),  i.e., 


P  = 


4RTsrour 


(50) 


To  test  the  predictions  of  this  part  of  the  model  apart  from  imperfections  in  the  predicted  value  of 
the  heat  of  vaporization,  I  use  the  experimental  value  of  Hvap  in  the  vapor-pressure  formulas.  With 
no  adjustment  parameters  at  all,  the  standard  deviation  of  the  predicted  from  experimental  values 
is  95%.  Evidently,  there  are  more  serious  shortcomings  in  the  vapor-pressure  model.  I 
experimented  with  several  empirical  modifications  of  the  model  and  got  interesting  results  by 
using  an  adjustable  parameter  to  scale  the  value  of  Hvap  in  computing  the  escape  velocity  of 
equation  49.  This  strategy  results  in  a  standard  deviation  of  about  36%,  a  much  improved 
accuracy  but  possibly  not  sufficient  for  use  in  the  multicomponent  evaporation  code.  The  best-fit 
scaling  factor  has  a  value  of  0.68,  i.e.,  the  calculation  is  significantly  improved  by  assuming  that 
only  68%  of  the  full  heat  of  vaporization  must  be  overcome  in  order  to  escape  the  liquid  surface. 
These  results  are  illustrated  in  Figure  22.  Noting  the  slight  downward  tendency  of  the  error  with 
increasing  Hvap  in  the  figure,  I  tried  using  a  two-parameter  fit  to  the  fraction  of  Hvap  used  to 
compute  ve.  This  improved  the  standard  deviation  slightly  to  30%.  The  model  may  well  be  further 
improved  and  placed  on  a  more  sound  theoretical  basis  by  doing  MD  studies  to  help  inform  the 
assumptions.  For  example,  perhaps  decreasing  the  minimum  escape  velocity  improves  the 
idealized  model  because  molecules  tend  to  equilibrate,  on  average,  at  a  value  of  potential  energy 
somewhat  above  that  in  the  bulk  by  means  of  collisions  in  the  interfacial  region  closest  to  the  bulk 
liquid.  A  great  advantage  of  this  bootstrapping  partnership  between  discrete  and  continuum 
descriptions  is  that,  by  using  the  same  model  potential  in  both,  the  physics  of  the  evaporation 
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Figure  22.  Accuracy  of  the  simple  vapor-pressure  model  using  61  L-J  parameter  sets  for  both  polar 
and  nonpolar  molecules. 

process  can  be  studied  and  built  into  the  continuum  model  apart  from  the  behavior  of  any  real 
substance.  Separate  studies  can  then  address  the  issue  of  the  best  potential  model  to  use.  It  goes 
without  saying  that  the  continuum  models  described  previously  can  be  implemented  with  any 
potential-energy  function  while  still  preserving  the  relatively  rapid  computational  qualities 
essential  to  incorporation  into  a  multicomponent  buming-rate  model. 

One  might  argue  that  empirically  based  engineering  correlations  could  do  considerably  better  at 
predicting  vapor  pressures.  In  fact,  I  applied  the  Reid  method  (44,  45)  based  on  the 
corresponding-states  principle  to  the  same  database  of  molecules  used  previously  and  found  that 
the  standard  deviation  was  only  4%.  However,  these  methods  are  apparently  much  less  successful 
when  applied  to  molecule  mixtures.  This  may  be  an  important  limitation  for  energetic-material 
combustion,  where  the  surface  is  multicomponent  even  for  the  simplest  case  of  ozone  as  explained 
in  section  3.3.2. 


43 


5.2  MD  Simulations  of  the  Condensed  Phase 

Another  fruitful  area  where  the  MD  approach  has  played  a  role  and  enjoys  bright  future  prospects 
is  in  the  condensed  phase.  An  example  of  past  use  is  in  determining  the  idealized  mass  densities 
(36)  of  putatively  pure  polymers  of  cellulose  mono-,  di-,  and  trinitrate.  These  values  are  needed  in 
the  semi-empirical  buming-rate  model  previously  discussed  for  nitrate-ester  propellants.  In 
addition,  the  method  was  also  applied  (36)  to  the  unreacting  solid  form  of  JA2  propellant  to  check 
the  accuracy  of  the  JA2  mass  density  computation  by  additive  molar  volumes,  the  method 
employed  by  the  CYCLOPS  (29,  30)  code.  Figure  23  shows  a  particular  relaxed  configuration  of 
JA2,  believed  to  be  the  first  published  “anatomically  correct”  molecular  view  of  a  real  propellant. 
The  CYCLOPS-computed  value  of  the  density  is  1.56  g/cm  ,  which  compares  very  well  with  the 
MD-computed  value  of  1.59  ±  0.02  g/cm  ,  which,  in  turn,  compares  very  well  with  experimental 
value  of  1.57  ±  0.01  g/cm3. 
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Figure  23.  One  view  of  a  relaxed  configuration  of  JA2  propellant  computed  by  a  MD 
simulation  (36)  consisting  of  2  chains  of  15  monomers  representing  NC,  16 
molecules  of  diethylene  glycol  dinitrate  (DEGDN),  and  10  molecules  of  NG. 
Oxygen  atoms  are  in  red,  nitrogen  in  blue,  carbon  in  grey,  and  hydrogen  in  white. 
A  number  of  the  component  molecules  can  be  identified  as  indicated.  This  is 
believed  to  be  the  first  computed  molecular  representation  of  a  real  propellant 
formulation. 
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Besides  the  purpose  just  described,  such  a  MD  simulation  could  be  used  to  determine  propellant 
heats  of  formation  that  properly  account  for  molecular  interactions  between  dissimilar  ingredients. 
Currently,  all  such  interactions  are  ignored  for  lack  of  computational  capability  and  the  expectation 
that  heats  of  solution  will  generally  be  small,  but  we  previously  found  (29,  30)  that  small  changes 
in  the  heats  of  formation  of  propellant  ingredients  can  have  surprisingly  large  effects  on  the 
burning  rate  at  low  and  intermediate  pressures.  This  effect  disappears  at  maximum  gun  pressures, 
and  therefore  will  not  affect  thermodynamic  equilibrium  calculations,  which  presently  must  make 
the  same  assumption  (that  heats  of  mixing  are  negligible).  The  reasons  for  this  unexpected 
sensitivity  of  the  burning  rate  to  ingredient  heats  of  formation  cannot  be  neatly  isolated  because  it 
is  a  highly  integrated  effect,  but  it  probably  arises  from  a  close  competition  between  energy  release 
by  reactions  and  the  locally  dissipative  processes  of  convection,  conduction,  and  molecular 
diffusion. 

Finally,  I  will  express  my  belief  that,  while  a  full  MD  description  of  the  gas  phase  will  never  be 
able  to  compete  with  the  continuum  description  in  terms  of  calculational  efficiency  and  accuracy 
(for  reasons  discussed  at  the  beginning  of  section  5),  neither  will  a  continuum  description 
ultimately  be  able  to  compete  with  a  discrete  molecular  description  of  the  condensed-phase 
processes.  It  may  be  a  long  time  before  this  promise  is  fulfilled,  but  if  a  complete  description  of 
propellant  combustion  is  to  be  realized,  there  appears  to  be  little  choice  but  to  pursue  the  MD 
approach  tenaciously.  In  my  discussions  with  molecular  dynamicists,  I  have  often  observed  an 
irrepressible  optimism  that  the  method  can  be  made  to  work  for  phase  changes,  for  subtle  transport 
effects  such  as  thermal  diffusion,  and  even  for  reactive  processes.  However,  examples  given  are 
usually  for  systems  that  are  idealized  in  the  extreme.  To  be  applicable  to  a  practical  propellant 
buming-rate  model,  these  simulations  will  have  to  be  made  to  work  for  very  general  systems,  even 
particularly  difficult  cases  involving  very  large  numbers  of  particles.  Perhaps  it  is  best  that  the 
dynamicists  not  fully  understand  the  difficulties  that  await  them! 


6.  Conclusions 


In  this  work,  I  have  made  a  reasoned  attempt  to  predict  the  most  promising  course  of  future 
research  aimed  at  predicting  the  burning  rates  of  solid  propellants  from  their  ingredients.  The  last 
15  years  has  seen  the  rise  to  dominance  of  models  that  treat  the  gas  phase  on  the  level  of 
elementary  reactions  with  full  multicomponent  transport.  This  explicit  recognition  of  chemical 
specificity  has  been  a  necessary  precursor  to  predictive  capability.  Energetic  materials  used  as 
propellants  are  designed  to  produce  gas  pressure  to  accomplish  work,  and  any  model  of  burning 
rate  must  include  descriptions  of  the  condensed  phase  and  the  gasification  mechanism.  However, 
attempts  to  treat  the  condensed-phase  and  surface  processes  at  the  same  level  of  rigor  as  the  gas 
phase  have  been  stymied  by  formidable  difficulties  in  both  experimental  and  theoretical 
approaches.  It  is  possible  that  experiments  may  never  be  devised  to  provide  the  kind  of  accuracy 
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and  detail  that  we  have  come  to  expect  in  the  gas  phase;  at  the  least,  it  will  probably  be  a  very  long 
time  in  coming.  Despite  this  sobering  possibility,  we  have  shown  that  some  degree  of 
predictability  may  be  possible  using  the  semi-empirical  approach  embodied  in  the  CYCLOPS 
code.  At  present,  this  predictability  appears  to  be  limited  to  members  of  well-studied  classes  of 
propellants,  such  as  double-base  and  RDX/HMX/binder  propellants.  However,  as  continuing 
research  encompasses  ever  wider  classes  of  ingredients,  generalizations  for  the  semi-empirical 
aspects  of  the  code  may  well  enable  wider  applicability.  For  example,  it  has  been  suggested  that  a 
single  pyrolysis  law  for  most  propellants  might  prove  sufficient  if  the  sensitivity  of  the  burning 
rate  to  this  law  is  low. 

In  my  opinion,  the  greatest  hope  for  treating  the  condensed-phase  and  surface  processes  in  a  full 
3-phase,  first-principles  model  lies  with  the  developing  field  of  reactive  MD.  Many  obstacles, 
both  known  and  unknown,  will  have  to  be  overcome  in  treating  condensed-phase  reactions  by  this 
approach.  It  is  not  clear  that  reactive  force  fields  with  sufficient  generality  can  be  developed.  It  is 
not  clear  that  methods  for  treating  the  many-body  interactions  can  be  developed  with  sufficient 
generality.  However,  there  are  ideas  in  the  community  about  how  to  proceed;  only  time  and  effort 
will  prove  their  value.  It  is  worth  remembering  that  burning  rate  is  a  highly  integrated 
macroscopic  consequence  of  almost  unfathomably  numerous  microscopic  processes.  Because  of 
this,  the  burning  rate  has  often  proved  to  be  extraordinarily  insensitive  to  the  underlying  processes. 
Thus,  one  may  be  justifiably  hopeful  that  even  an  imperfect  description  of  the  detailed  processes 
may  lead  to  predictions  of  macroscopic  phenomena  that  will  provide  insights  and  guidance  of  a 
practical  nature  in  the  development  of  new  and  optimized  propellant  formulations. 
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Appendix.  Propellant  Ingredients  and  Formulations 


AMMO  poly  3-azidomethyl-3-methyl  oxetane 
BAMO  poly  3,3-(bis)azidomethyl  oxetane 

CBIH  copolymer  of  butadiene  and  isoprene  with  hydroxyl  terminated  groups 

CL20  2.4. 6. 8. 1 0, 1 2-hexanitrohexaazaisowurtzitane,  rectangular  crystal  size  5x3x3  pm 

DEGDN  diethylene  glycol  dinitrate 

GAP1U  glycidyl  azide -poly urethane  copolymer 

GAP2  glycidyl  azide  polymer  (molecular  mass  of  2000) 

HMX  cyclotetramethylenetetranitramine 

PU  polyurethane 

PUNE  mixture  of  PU  and  NE 


NC  nitrocellulose  (cellulose  nitrate) 

NCI  cellulose  mononitrate 


NC2  cellulose  dinitrate 


NC3  cellulose  trinitrate 


NE  mixture  of  two  nitroesters  (dinitratdietileneglycole  and  dinitrattrietileneglicole) 

NG  nitroglycerin  (glycerin  trinitrate) 

RDX  cyclotrimethylenetrinitramine 

THF  tetrahydrofuran 


Nominal  Propellant  Formulations  Assumed  by  CYCLOPS  Code: 
JA2  60%  NC  (13.1  %N),  25%  DEGDN,  15%  NG 
M9  59.1%  NC  (13.25  %N),  40.9%  NG 
M10  100%  NC  (13.15  %N) 
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Intentionally  left  blank. 
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ABERDEEN  PROVING  GROUND 
DIR  USARL 

AMSRL  Cl  LP  (BLDG  305) 
AMSRL  Cl  OK  TP  (BLDG  4600) 


1  COMMANDING  GENERAL 
US  ARMY  MATERIEL  CMD 
AMCRDA  TF 
5001  EISENHOWER  AVE 
ALEXANDRIA  VA  22333-0001 

1  INST  FOR  ADVNCD  TCHNLGY 
THE  UNIV  OF  TEXAS  AT  AUSTIN 
3925  W  BRAKER  LN  STE  400 
AUSTIN  TX  78759-5316 

1  US  MILITARY  ACADEMY 

MATH  SCI  CTR  EXCELLENCE 
MADN  MATH 
THAYER  HALL 
WEST  POINT  NY  10996-1786 

1  DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRL  D 
DR  D  SMITH 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1197 

1  DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRL  CS  IS  R 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1197 

3  DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRL  Cl  OK  TL 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1197 

3  DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRL  CS  IS  T 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1197 


2DEFENSE  TECHNICAL 

INFORMATION  CENTER 

DTIC  OCA  2 

8725  JOHN  J  KINGMAN  RD 
STE  0944 

FT  BEL  VOIR  VA  22060-6218 
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NO.  OF 

COPIES  ORGANIZATION 

ABERDEEN  PROVING  GROUND 


28  DIR  USARL 

AMSRL  WM  BD 
W  R  ANDERSON 
R  A  BEYER 
A  L  BRANT 
S  W  BUNTE 
C  F  CHABALOWSKI 
L  M  CHANG 
T  P  COFFEE 
J  COLBURN 
P  J  CONROY 
R  A  FIFER 
B  E  FORCH 
B  E  HOMAN 
S  L  HOWARD 
P  J  KASTE 
AJKOTLAR 
C  LEVERITT 
K  L  MCNESBY 
M  MCQUAID 
M  S  MILLER 
A  W  MIZIOLEK 
J  B  MORRIS 
J  A  NEWBERRY 
M  J  NUSCA 

R  A  PESCE-RODRIGUEZ 
G  P  REEVES 
B  M  RICE 
R  C  SAUSA 
A  W  WILLIAMS 
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NO.  OF 

COPIES  ORGANIZATION 

1  M  MILLER 

1124  COWPENS  AVE 
TOWSON  MD  21286 
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